//
//  LIBRERIA DI ALGEBRA LINEARE - VERSIONE 20080828 
//
var versione_alin= 20080828;
//
//<![CDATA[
//
//   L'equivalente in Fortran 90 del type(ARMA), a parte la stringa commento,
//   è definito in questo modo:
//
//     type,public:: arma
//         integer::error,ordine,ld
//         real,dimension(:),pointer::a  ! vettore di (ld+1)*(ordine+1) elementi. 
//         real::scalare
//         integer::clip_sr,clip_srold,clip_sc,clip_scold
//         integer::clip_er,clip_erold,clip_ec,clip_ecold
//         integer::ndim
//     end type 
//
//   Ogni elemento di type(ARMA) va opportunamente inizializzato in Fortran 90
//   seguendo le regole indicate nella corrispondente function ARMA in
//   javascript.
//
function ARMA(ld,n){
    var j,k,i;
    // Questo tipo derivato serve per definire sia vettori che matrici.
    //
    // Commento a questo dato:
    this.commento="<br/>"
    // Risultato dell'ultima operazione (0 == riuscita bene).
    this.error=0;
    // Ordine della matrice rettangolare, conservato per sicurezza.
    // Se la matrice e' un vettore deve valere 0.
    this.ordine=parseInt(n);
    // Lunghezza di ogni colonna meno 1 (ossia quante righe stanno
    // in una colonna). 
    this.ld=parseInt(ld);
    if(this.ld < this.ordine) this.ld=this.ordine;
    // Dati della matrice e della soluzione.
    this.a= new Array();
    //
    // In 0 ci sta l'ordine del sistema ossia n
    // Tra 1 e n ci sta, a volte, il vettore soluzione.
    // Ogni colonna e' lunga ld >=n e l'elemento 0
    // di ogni colonna serve per gestire i dati ed e' 
    // un intero mentre gli altri sono float.
    //
    i=-1;
    for (j=0;j<=this.ordine;j++){
        i++;
        this.a[i]=j;
        for(k=1;k<=this.ld;k++){
            i++;
            this.a[i]=Math.PI;       
            }
        }
    this.a[0]=this.ordine;
    // Eventuale risultato scalare di un calcolo algebrico.
    this.scalare=Math.PI;
    // Dati per definire la sottomatrice di lavoro.
    // Primo elemento della sottomatrice (start riga,colonna).
    this.clip_sr=1;
    this.clip_srold=this.clip_sr;
    this.clip_sc=1;
    this.clip_scold=this.clip_sc;
    // Ultimo elemento della sottomatrice (end riga,colonna).
    // this.clip_er=this.ordine;
    if(this.ordine>0) { this.clip_er=this.ordine; }
    else { this.clip_er=this.ld;}
    this.clip_erold=this.clip_er;
    this.clip_ec=this.ordine;
    this.clip_ecold=this.clip_ec;
    // Limite del vettore ed ordine delle matrici di servizio
    this.ndim=this.ordine;
    if(this.ordine==0) this.ndim=this.ld;
    return this;
    }
//_________________________________________________________________
//
//  Funzioni che usano le sottomatrici delimitate dal CLIP.
//  I nomi iniziano tutti con clip.
//
function clip(a,sr,sc,er,ec){
    //
    // Definisce i limiti di clip e verifica che siano corretti.
    // I valori vecchi sono recuperabili con la clip_old.
    //
    // a.clip_sr == start riga ( da 1 in poi )
    // a.clip_sc == start colonna ( da 0 in poi )
    // a.clip_er == end riga ( da 1 in poi )
    // a.clip_ec == end colonna (da 0 in poi )
    // in colonna 0 ci sta il vettore di servizio.
    //
    a.clip_srold=a.clip_sr;
    a.clip_scold=a.clip_sc;
    a.clip_erold=a.clip_er;
    a.clip_ecold=a.clip_ec;
    a.clip_sr=parseInt(sr);
    a.clip_sc=parseInt(sc);
    a.clip_er=parseInt(er);
    a.clip_ec=parseInt(ec);
    if(a.clip_sr<1)a.clip_sr=1;
    // Consentito operare anche sulla colonna 0
    if(a.clip_sc<0)a.clip_sc=1;
    if(a.clip_er<1)a.clip_er=1;
    if(a.clip_ec<0)a.clip_ec=1;
    if(a.clip_sr>a.ld)a.clip_sr=a.ld;
    if(a.clip_sc>a.ordine)a.clip_sc=a.ordine;
    if(a.clip_er>a.ld)a.clip_er=a.ld;
    if(a.clip_ec>a.ordine)a.clip_ec=a.ordine;
    if(a.clip_er<a.clip_sr) a.clip_er=a.clip_sr;
    if(a.clip_ec<a.clip_sc) a.clip_ec=a.clip_sc;
    }
//
function clip_coerenza(a){
    //
    //  I dati di clip devono individuare sempre almeno
    //  un elemento realmente contenuto nella matrice.
    //
    if(a.clip_sr<1)a.clip_sr=1;
    if(a.clip_sc<0)a.clip_sc=1;
    if(a.clip_er<1)a.clip_er=1;
    if(a.clip_ec<0)a.clip_ec=1;
    if(a.clip_sr>a.ld)a.clip_sr=a.ld; // vietato superare a.ld ossia la leading dimension
    if(a.clip_sc>a.ordine)a.clip_sc=a.ordine; // vietato superare l'ordine
    if(a.clip_er>a.ld)a.clip_er=a.ld;
    if(a.clip_ec>a.ordine)a.clip_ec=a.ordine;
    if(a.clip_er<a.clip_sr) a.clip_er=a.clip_sr;
    if(a.clip_ec<a.clip_sc) a.clip_ec=a.clip_sc;
    }
//
function clip_commento(a){
    // 
    //  Aggiorna il commento se è stato indicato e fornisce come 
    //  risultato il commento ( esistente prima dell'eventuale aggiornamento).
    //
    var lista=clip_commento.arguments;
    var str=a.commento;
    if(lista.length>1) a.commento=lista[1];
    return str;
    }
//
function clip_scalare(a){
    //
    // Aggiorna il valore dello scalare interno e fornisce come risultato
    // il vecchio valore (che puo' essere non aggiornato se il nuovo non
    // e' stato indicato com secondo argomento della lista di chiamata )
    //
    var lista=clip_commento.arguments;
    var sc=a.scalare;
    if(lista.length>1) a.scalare=lista[1];
    return sc;
    }
//
function clip_get(a){
    // 
    // Restituisce un vettore con i dati dell'oggetto ARMA
    //
    return [ "Commento:"+a.commento+
             a.ld+"=leading dimension ",
             a.ordine+"=ordine ",
             a.ndim+"=ordine matrici di servizio ",
             a.scalare+"=valore scalare memorizzato ",
             a.clip_sr+"=riga iniziale clip ",
             a.clip_er+"=riga finale clip ",
             a.clip_sc+"=colonna iniziale clip ",
             a.clip_ec+"=colonna finale clip ",
             a.error+"=indice di errore" ];
    }
//
function clip_copymat(a,b){
    //
    // Copiatura tra sottomatrici definite dalla clippatura.
    //
    var ma=a.ld+1;
    var mb=b.ld+1;
    var deltar,deltac,bdifr,bdifc;
    var j,k,js,je,ks,ke;
    deltar=a.clip_er-a.clip_sr;
    deltac=a.clip_ec-a.clip_sc;
    bdifr=b.clip_sr-a.clip_sr;
    bdifc=b.clip_sc-a.clip_sc;
    j=b.clip_er-b.clip_sr;
    k=b.clip_ec-b.clip_sc;
    if(j<deltar)deltar=j;
    if(k<deltac)deltac=k;
    js=a.clip_sc;
    je=a.clip_sc+deltac;
    ks=a.clip_sr;
    ke=a.clip_sr+deltar;
    //
    for(j=js;j<=je;j++){
        for(k=ks;k<=ke;k++){
            a.a[k+j*ma]=b.a[k+bdifr+(j+bdifc)*mb];
            }
        }
    }
//
function clip_diade(a,b,c){
    //
    // La prima colonna della sottomatrice b moltiplicata per
    // la prima riga della sottomatrice c formano una diade
    // aggiunta alla sottomatrice a.
    //
    var j,k,bcost,ccost,deltar,deltac,js,je,ks,ke,bdifr,cdifc
    ma=a.ld+1;
    mb=b.ld+1;
    mc=c.ld+1;
    j=a.clip_ec-a.clip_sc;
    deltac=c.clip_ec-c.clip_sc;
    if(deltac>j) deltac=j;
    j=a.clip_er-a.clip_sr;
    deltar=b.clip_er-b.clip_sr;
    if(deltar>j) deltar=j;
    js=a.clip_sc;
    je=js+deltac;
    ks=a.clip_sr;
    ke=ks+deltar;
    bdifr=b.clip_sr-a.clip_sr;
    cdifc=c.clip_sc-a.clip_sc;
    //
    bcost=bdifr+b.clip_sc*mb;
    ccost=c.clip_sr+cdifc*mc;
    for(j=js;j<=je;j++){
        for(k=ks;k<=ke;k++){
            a.a[k+j*ma]+=b.a[k+bcost]*c.a[j*mc+ccost];
            }
        }
    }
//
function clip_diagomat(a,diago,nondiago){
    //
    // La diagonale della sottomatrice vale il numero diago, il resto
    // vale il numero nondiago.
    //
    var i,j,k;
    for(j=a.clip_sc;j<=a.clip_ec;j++) {
        i=j*(a.ld+1);
        for(k=a.clip_sr;k<=a.clip_er;k++){
            a.a[i+k]=nondiago;
            }
        a.a[i+j+a.clip_sr-a.clip_sc]=diago;
        }
    return 0;
    }
//
// ATTENZIONE: impostare xfatondo ad un valore nullo se non
// si vogliono arrotondamenti !
//
//
var xumfatondo=1.e-8
var xfatondo=parseInt(1/xumfatondo);
function fatondo(a){ // arrotonda un numero reale
   var b,c;
   if(a>=0){
          b=a+xumfatondo/2;
          c=parseInt(b*xfatondo)*xumfatondo;
       }
   else {
          b=-a+xumfatondo/2;
          c=-parseInt(b*xfatondo)*xumfatondo;
       }
   return c;
   }
//
var clip_formato=12; // Usare almeno 5 (cinque cifre utili)
function clip_format(aa){ // arrotonda un numero reale
   var a,b,c,n,s;
   var s=aa+"                         ";
   n=s.indexOf("e");
   if(n<0) a= parseFloat(s.substr(0,clip_formato));
   else a=aa;
//   alert("s= "+s+" n= "+n+"  a= "+a);
   return a;
   }
//
function clip_getmat(v,sinfo){
    //
    // Preleva valori dalla sottomatrice "clippata" entro v, tramite una 
    // stringa sinfo contenente dati float separati tra loro dal 
    // carattere identico all' ultimo carattere della stringa.
    // L'inizio della stringa deve essere un intero che 
    // specifica quanti dati vengono prelevati al massimo.
    // Per esempio sinfo potrebbe essere "8|" ossia richiedere 8
    // dati e chiedere che siano separati dal carattere "|".
    // 
    // Per inserire dati in v usare la clip_setmat(...).
    //
    var j,k,nn,ne,ni,nf,nj,mo,me,mm,form,separatore;
    var s="";
    if (sinfo.length > 0) {
        separatore=sinfo.substr(sinfo.length-1);
        }
    else {
        v.error=-1;
        return s;
        }
    quanti=parseInt(sinfo);
    if(isNaN(quanti)){
        v.error=-777;
        return s;
        }
    nn=v.clip_sr;
    ne=v.clip_er;
    mo=v.clip_sc;
    me=v.clip_ec;
    ni=0;
    for (j=mo; j<=me; j++){
        nj=nn+j*(v.ld+1);
        mm=nj+ne-nn;
        for (k=nj;k<=mm;k++){
            ni++;
            if(ni>quanti) break;
            form=v.a[k]+" ";
            s=s+form+separatore;
            }
        }   
    v.error=0;
    //
    //  Ritorna la stringa 
    //
    return s;       
    }
//
function clip_inversa(a,usemat){
    //
    // Fa l'inversa sostituendola alla matrice originaria.
    // Come ordine usa la differenza delle colonne piu' uno.
    // ATTENZIONE:
    // Usa una matrice esterna di servizio (usemat) che deve essere di 
    // dimensioni adeguate.
    //
    var j,base,riesco,ordine=a.clip_ec-a.clip_sc+1
    if(ordine > a.clip_ec-a.clip_sc+1) {
        a.error=-101;
        return a.error;
        } 
    if(isNaN(usemat)){
        if(ordine > usemat.ordine) {
            alert("Inversa non fattibile\n"+
              "L'ordine della matrice di servizio vale "+usemat.ordine);
            a.error=-102;
            return a.error;
            }
         }
    else {
       alert("usemat risulta essere un numero e non un oggetto di tipo ARMA"+
             "\nInizializzare usemat !");
         a.error=-1000;
         return a.error;
         }
    base=a.clip_sc-1;
    usemat.clip_sc=1;
    usemat.clip_sr=1;
    usemat.clip_ec=ordine;
    usemat.clip_er=ordine;
    usemat.ndim=ordine;
    clip_copymat(usemat,a);
    usemat.clip_sc=0;
    usemat.clip_ec=0;
    riesco=lu_fa(usemat);
    if(riesco!=0) {
        alert("fattorizzazione fallita in clip_inversa: "+riesco);
        a.error=riesco;
        return a.error;
        }
    clip_save(a);
    for(j=1; j<=ordine;j++){
        zerovet(usemat);
        usemat.a[j]=1.0;
        lu_so(usemat);
        a.clip_sc=base+j;
        a.clip_ec=base+j;
        clip_copymat(a,usemat);
        }
    clip_old(a);
    clip_old(usemat);
    }
//
function clip_matmul(a,b,c){
    //
    // Prodotto matriciale tra le sottomatrici b e c che
    // viene aggiunto alla sottomatrice a.
    //
    var j,k,bcost,ccost,deltar,deltac,js,je,ks,ke,bdifr,cdifc;
    var p,passi;
    ma=a.ld+1;
    mb=b.ld+1;
    mc=c.ld+1;
    j=a.clip_ec-a.clip_sc;
    deltac=c.clip_ec-c.clip_sc;
    if(deltac>j) deltac=j;
    j=a.clip_er-a.clip_sr;
    deltar=b.clip_er-b.clip_sr;
    if(deltar>j) deltar=j;
    js=a.clip_sc;
    je=js+deltac;
    ks=a.clip_sr;
    ke=ks+deltar;
    bdifr=b.clip_sr-a.clip_sr;
    cdifc=c.clip_sc-a.clip_sc;
    passi=b.clip_ec-b.clip_sc;
    j=c.clip_er-c.clip_sr;
    if(passi >j) passi=j;
    //
    for(p=0;p<=passi;p++){
        bcost=bdifr+(b.clip_sc+p)*mb;
        ccost=(c.clip_sr+p)+cdifc*mc;
        for(j=js;j<=je;j++){
            for(k=ks;k<=ke;k++){
                a.a[k+j*ma]+=b.a[k+bcost]*c.a[j*mc+ccost];
                }
            }
        }
    }
//
function clip_matpluseq(a,scala,b){
    //
    // Alla sottomatrice a aggiunge la sottomatrice b scalata.
    //
    var ma=a.ld+1;
    var mb=b.ld+1;
    var deltar,deltac,bdifr,bdifc;
    var j,k,js,je,ks,ke;
    deltar=a.clip_er-a.clip_sr;
    deltac=a.clip_ec-a.clip_sc;
    bdifr=b.clip_sr-a.clip_sr;
    bdifc=b.clip_sc-a.clip_sc;
    j=b.clip_er-b.clip_sr;
    k=b.clip_ec-b.clip_sc;
    if(j<deltar)deltar=j;
    if(k<deltac)deltac=k;
    js=a.clip_sc;
    je=a.clip_sc+deltac;
    ks=a.clip_sr;
    ke=a.clip_sr+deltar;
    //
    for(j=js;j<=je;j++){
        for(k=ks;k<=ke;k++){
            a.a[k+j*ma]+=scala*b.a[k+bdifr+(j+bdifc)*mb];
            }
        }
    }
//
function clip_matscala(a,scala){
    //
    //   Effettua la scalatura della matrice a (la sola sottomatrice clippata)
    //
    var ma=a.ld+1;
    var deltar,deltac;
    var j,k,js,je,ks,ke;
    deltar=a.clip_er-a.clip_sr;
    deltac=a.clip_ec-a.clip_sc;
    js=a.clip_sc;
    je=a.clip_sc+deltac;
    ks=a.clip_sr;
    ke=a.clip_sr+deltar;
    //
    for(j=js;j<=je;j++){
        for(k=ks;k<=ke;k++){
            a.a[k+j*ma]*=scala;
            }
        }
    }
//
function clip_old(a){
    //
    // Ristabilisce i vecchi limiti di clip
    //
    a.clip_sr=a.clip_srold;
    a.clip_sc=a.clip_scold;
    a.clip_er=a.clip_erold;
    a.clip_ec=a.clip_ecold;
    }
//
function clip_save(a){
    //
    // Salva i limiti di clip senza aggiornarli. Con la clip
    // oltre al salvataggio dei vecchi avviene l'aggiornamento.
    //
    a.clip_srold=a.clip_sr; // start riga
    a.clip_scold=a.clip_sc; // start colonna
    a.clip_erold=a.clip_er; // end riga
    a.clip_ecold=a.clip_ec; // end colonna
    }
//
function clip_setmat(v,s){
    //
    // Assegna valori alla sottomatrice "clippata" tramite una 
    // stringa contenente dati float separati tra loro dal 
    // carattere identico all' ultimo carattere della stringa.
    //
    // Per estrarre dati da v usare la clip_getmat(...).
    // 
    // Attenzione: i dati vanno dati per colonna ossia
    // prima la prima colonna, poi la seconda etc...
    // Se mancano dati per riempire la matrice restano invariati
    // i vecchi dati.
    // I dati in eccesso sulla stringa vengono ignorati.
    //
    var j,k,nn,ne,ni,nf,nj,mo,me,mm,separatore;
    if (s.length > 0) {
        separatore=s.substr(s.length-1);
        }
    else {
        v.error=-10;
        return v.error;
        }
    nn=v.clip_sr;
    ne=v.clip_er;
    mo=v.clip_sc;
    me=v.clip_ec;
    //alert("Separatore "+separatore);
    ni=0;
    for (j=mo; j<=me; j++){
        nj=nn+j*(v.ld+1);
        mm=nj+ne-nn;
        for (k=nj;k<=mm;k++){
            nf=s.indexOf(separatore,ni);
            if(nf<0) break;
            v.a[k]=parseFloat(s.slice(ni,nf));
            ni=nf+1;
            }
        }   
    v.error=0;
    //
    //  Ritorna il punto da cui riprendere per estrarre 
    //  altri dati dalla stringa.
    //
    return ni;       
    }
//
function clip_trasposta(a,b){
    //
    // Copiatura tra sottomatrici. In a viene copiata la trasposta di b.
    //
    var ma=a.ld+1;
    var mb=b.ld+1;
    var deltar,deltac,bdifr,bdifc;
    var j,k,js,je,ks,ke;
    deltar=a.clip_er-a.clip_sr;
    deltac=a.clip_ec-a.clip_sc;
    bdifr=b.clip_sr-a.clip_sc;
    bdifc=b.clip_sc-a.clip_sr;
    j=b.clip_er-b.clip_sr;
    k=b.clip_ec-b.clip_sc;
    if(k<deltar)deltar=k;
    if(j<deltac)deltac=j;
    js=a.clip_sc;
    je=a.clip_sc+deltac;
    ks=a.clip_sr;
    ke=a.clip_sr+deltar;
    //
    for(j=js;j<=je;j++){
        for(k=ks;k<=ke;k++){
            a.a[k+j*ma]=b.a[j+bdifr+(k+bdifc)*mb];
            }
        }
    }
//
var clip_maxcol=15;  // Non stampa piu' di queste colonne.
function clip_vedomat(a){
    //
    // Fa la parte interna di una tabella ossia usa solo i tag <tr> e <td>
    // Per stampare il risultato aggiungere "<table>" e concludere
    // con </table>.
    //
    var j,k,ma,js,je,ks,ke,ss="";
    ma=a.ld+1;
    js=a.clip_sc;
    je=a.clip_ec;
    ks=a.clip_sr;
    ke=a.clip_er;
    if(je-js > clip_maxcol) je=js+clip_maxcol;
    for(k=ks; k<=ke; k++){
        ss=ss+"<tr>";
        for(j=js; j<=je; j++){
            // ss=ss+"<td><input type=text value='"+a.a[k+ma*j]+"' size=5></td>";
            ss=ss+"<td>"+clip_format(fatondo(a.a[k+ma*j]))+"</td>"
            }
        ss=ss+"</tr>";
        }
    return ss;
    }
//
function clip_vedovet(a){
    //
    // Fa la parte interna di una tabella ossia usa solo i tag <tr> e <td>
    // Visualizza i dati come vettore colonna.
    // Per stampare il risultato aggiungere "<table>" e concludere
    // con </table>.
    //
    // Visualizza la colonna 0 di a ossia il vettore di servizio.
    //
    var k,ks,ke,ss="";
    ks=a.clip_sr;
    ke=a.clip_er;
    for(k=ks; k<=ke; k++){
        ss=ss+"<tr><td>"+k+")</td><td>"+clip_format(fatondo(a.a[k]))+"</td></tr>"
        }
    return ss;
    }
//
//________________________________________________________________________
//
//  Funzioni che NON usano i dati di clip e quindi operano usando
//  il componente matrice.ndim.
//
//  Queste funzioni vengono tratte da libri di matematica e dunque
//  rappresentano la parte realmente "matematica" di questa libreria.
//
function copyvet(va,vb){
    // NON usa i limiti di clip ma va da 1 a va.ndim o vb.ndim
    // Copia la colonna 0 di vb nella colonna 0 di va ossia
    // trasferisce il vettore di servizio dell'oggetto ARMA vb 
    // nel vettore di servizio dell'oggetto ARMA va.
    var j,mm=vb.ndim;
    if(mm>va.ndim)mm=va.ndim;
    for(j=1;j<=mm;j++){
        va.a[j]=vb.a[j];
        }
    }
//
//  Uno dei modi per inserire elementi in un oggetto ARMA è quello di usare
//  il suo componente a come un normale vettore e poi posizionare il vettore
//  di servizio sulla riga o colonna che serve.
//  Ecco le function da utilizzare...
//
function caricavet(a,vest){
    // vest è un normale vettore esterno che viene usato in due modi.
    // Se l'elemento 0 di vest ossia v[0] e' non negativo ( 0 o positivo) allora gli 
    // altri dati vengono trasferiti nel vettore di servizio di a.
    // Se l'elemento 0 di vest e' negativo allora vengono prelevati
    // gli elementi dal vettore di servizio di a ( che deve essere un
    // oggetto ARMA ) e messi in vest da 1 al minore tra vest.length-1
    // e a.ndim
    var j,mm;
    mm=vest.length;
    if(isNaN(mm)){
         alert("In caricavet : errore: secondo argomento non vettore");
         return -1;
         }
    mm-=1;
    if(mm>a.ndim) mm=a.ndim;
    if(1>mm){
         alert("In caricavet : nessuna operazione. Vettori troppo corti");
         return -2;
         }
    if(0>vest[0]){
        for(j=1;j<=mm;j++){
            vest[j]=a.a[j];
            }
         }
    else {
        for(j=1;j<=mm;j++){
            a.a[j]=vest[j];
            }
         }
    return 0;
    }
//
function copiariga(a,b,kk){
    // Se k e' non negativo ( ossia 0 o positivo) trascrive la colonna 0 di b 
    // nella riga k di a.
    // Se k e' negativo fa il contrario ossia trascrive la riga abs(k) di a 
    // nella colonna 0 di b.
    // Il trasferimento parte da 1 e prende come limite il minore 
    // tra a.ndim e b.ndim.
    var j,k,passo;
    var mm=a.ndim;
    if(mm>b.ndim) mm=b.ndim;
    if(0>kk){
         k=-kk;
         passo=a.ld+1;
         for(j=1;j<=mm;j++){
             k+=passo;
             b.a[j]=a.a[k];
             }
         }
    else {
         k=kk;
         passo=b.ld+1;
         for(j=1;j<=mm;j++){
             k+=passo;
             a.a[j]=b.a[k];
             }
         }
    }
//
function copiacolonna(a,b,kk){
    // Se k e' non negativo ( ossia 0 o positivo) trascrive la colonna 0 di b 
    // nella colonna k di a.
    // Se k e' negativo fa il contrario ossia trascrive la colonna abs(k) di a 
    // nella colonna 0 di b.
    // Il trasferimento parte da 1 e prende come limite il minore 
    // tra a.ndim e b.ndim.
    var j,k,passo;
    var mm=a.ndim;
    if(mm>b.ndim) mm=b.ndim;
    if(0>kk){
         passo=a.ld+1;
         k=-kk*passo;
         for(j=1;j<=mm;j++){
             b.a[j]=a.a[j+k];
             }
         }
    else {
         passo=b.ld+1;
         k=kk*passo;
         for(j=1;j<=mm;j++){
             a.a[j]=b.a[j+k];
             }
         }
    }
//
// I dati di un oggetto ARMA possono essere presi anche da ( o salvati in )
// un vettore di vettori ossia da un vettore trattato come vettore colonna, 
// i cui elementi sono considerati vettori riga.
// Se i singoli vettori riga sono di lunghezza diversa vengono usati
// al meglio. La numerazione di un vettore di vettori parte dallo
// zero mentre la matrice dell'oggetto ARMA e' numerata a partire da 1.
//
function davetdivet(a,vdv){
    // Riempie la matrice a prendendo i dati presenti nel vdv ossia
    // nel vettore di vettori.
    var j,k,kk,mm,nn,passo;
    passo=a.ld+1;
    mm=a.ld;
    if(mm>vdv.length) mm=vdv.length;
    for(j=1;j<=mm;j++){
        nn=a.ndim;
        if(nn>vdv[j-1].length) nn=vdv[j-1].length;
        kk=j;
        for(k=1;k<=nn;k++){
             kk+=passo;
             a.a[kk]=vdv[j-1][k-1];
             }
        }
    }
//
function invetdivet(a,vdv){
    // Riempie il  vdv ossia il vettore di vettori prendendo i dati dalla
    // matrice contenuta in a.
    var j,k,kk,mm,nn,passo;
    passo=a.ld+1;
    mm=a.ld;
    if(mm>vdv.length) mm=vdv.length;
    for(j=1;j<=mm;j++){
        nn=a.ndim;
        if(nn>vdv[j-1].length) nn=vdv[j-1].length;
        kk=j;
        for(k=1;k<=nn;k++){
             kk+=passo;
             vdv[j-1][k-1]=a.a[kk];
             }
        }
    }
//
function euclidea(v){
    // NON usa i limiti di clip. Norma euclidea del vettore v.
    var k;
    v.scalare=0.0;
    for(k=1;k<=v.nmax;k++){
        v.scalare+=v.a[k]*v.a[k];
        }
    v.scalare=Math.sqrt(v.scalare);
    v.error=0;
    return v.scalare;
    }
//
function getvet(v,s,primo,ultimo){
    // Produce una stringa contenente i valori del vettore di servizio
    // dell'oggetto v ( di tipo ARMA) dall'elemento primo all'elemento ultimo.
    // Ricordarsi che in s va messa una stringa il cui ultimo carattere
    // e' il separatore tra i dati e che deve iniziare con un numero
    // che specifica il massimo numero di dati da estrarre. Questo limite
    // prevale su quello dato da (ultimo-primo+1).
    // NON usa i limiti di clip
    var ss;
    clip(v,primo,0,ultimo,0)
    ss=clip_getmat(v,s);
    clip_old(v);
    return ss; // Attenzione: il risultato e' una stringa.
    }
//
function hessenberg(a){
    //  Attenzione: la trasformazione di  Hessemberg produce 
    //  una matrice in cui tutti gli elementi della triangolare
    //  inferiore, ad eccezione della prima sottodiagonale 
    //  sono nulli.
    //  Questa function non azzera gli elementi ma li usa
    //  per memorizzare i vettori della trasformazione.
    //  La trasformazione di Hessemberg conserva gli autovalori
    //  della matrice per cui il risultato puo' essere 
    //  passato alla function hqr che trova gli autovalori
    //  di una matrice in forma di Hessemberg (distruggendo
    //  i dati di servizio prodotti da questa function)
    //
    var i,j,m,x,y;
    var ld=a.ld+1;
    var n=a.ndim;
    for (m=2;m<n;m++){
        x=0.0;
        i=m;
        for(j=m;j<=n;j++){
            if(Math.abs(a.a[j+(m-1)*ld]) > Math.abs(x) ){
                 x=a.a[j+(m-1)*ld];
                 i=j;
                 }
            }
        if(i!=m){
            for(j=m-1;j<=n;j++){
                y=a.a[i+j*ld];
                a.a[i+j*ld]=a.a[m+j*ld];
                a.a[m+j*ld]=y;
                }
            for(j=1;j<=n;j++){
                y=a.a[j+i*ld];
                a.a[j+i*ld]=a.a[j+m*ld];
                a.a[j+m*ld]=y;
                }
            }
        if(Math.abs(x)>0.0) {
            for(i=m+1;i<=n;i++){
                y=a.a[i+(m-1)*ld];
                if(Math.abs(y)>0.0) {
                    y=y/x;
                    a.a[i+(m-1)*ld]=y;
                    for(j=m;j<=n;j++){
                        a.a[i+j*ld]+= -y*a.a[m+j*ld];
                        }
                    for(j=1;j<=n;j++){
                        a.a[j+m*ld]+=y*a.a[j+i*ld];
                        }
                    }
                }
            }
        }
    }
//
var hqr_itsmax=30;
function hqr(a){
    //
    //  L'algoritmo QR di J.C.F.Francis e V.N.Kublanovskaya.
    //  Gli autovalori di a sono memorizzati in a.a[1]...a.a[a.ndim]
    //  per la parte reale e in a.a[ld]...a.a[ld*a.ndim] per la 
    //  parte immaginaria dove ld=a.ld+1.
    //  E' prevista una apposita funzione ausiliaria per prelevare il 
    //  vettore della parte immaginaria ricopiandolo nel vettore voluto.
    //  Il risultato e' il numero di iterazioni fatte. Se e' negativo
    //  il calcolo e' stato interrotto per raggiunto limite di 
    //  iterazioni. 
    //
    //  Da Numerical Recipes in Fortran, II edition, Cap. 11 pag 484  
    //  "The QR Algorithm for Real Hessenberg Matrices" 
    //
    var ld,n,i,its,j,k,l,m,nn;
    var anorm,p,q,r,s,t,u,v,w,x,y,z;
    ld=a.ld+1;
    n=a.ndim;
    anorm=Math.abs(a.a[1+ld]);
    for (i=2; i<=n; i++) {
        for(j=i-1; j<=n;j++){
            anorm += Math.abs(a.a[i+j*ld]);
            }
        }
    // alert("Entra in hqr");
    nn=n;
    t=0.0;
    for(eigen=0;eigen<n;eigen++) {
        if(nn<1) break;
        // ciclo its
        // alert("autovalore: "+nn);
        for(its=0; its<4*hqr_itsmax;its++) {
            l=1;
            for(i=nn;i>1;i--){
                s=Math.abs(a.a[(i-1)+ld*(i-1)])+
                  Math.abs(a.a[i+ld*i])
                if(Math.abs(s)<=0.0) s=anorm;
                if((Math.abs(a.a[i+ld*(i-1)])+s) <= s) {
                    l=i; 
                    break;
                    }
                }
            // 3
            x=a.a[nn+ld*nn];
            if(l == nn) {
                // un solo autovalore
                a.a[nn]=x+t;
                a.a[nn*ld]=0.0;
                nn--;
                break; // esce dalle its
                }
            else {
                y=a.a[nn-1+ld*(nn-1)];
                w=a.a[nn+ld*(nn-1)]*a.a[nn-1+ld*nn];
                if(l==nn-1) {
                    // una coppia di autovalori
                    p=0.5*(y-x);
                    q=p*p+w;
                    z=Math.sqrt(Math.abs(q));
                    x=x+t;
                    if(q >=0.0) {
                        // Due autov. reali
                        // z=p+sign(z,p)       ATTENZIONE QUI.... z inizialmente positivo
                        if(p <0.0) z=p-z;
                        else z=p+z;
                        a.a[nn]=x+z;
                        a.a[nn-1]=a.a[nn];
                        if(Math.abs(z)>0.0) a.a[nn]=x-w/z;
                        a.a[ld*nn]=0.0;
                        a.a[ld*(nn-1)]=0.0;
                        }
                    else {
                        // Due autov. complessi coniug.
                        a.a[nn]=x+p;
                        a.a[nn-1]=a.a[nn];
                        a.a[ld*nn]=z;
                        a.a[ld*(nn-1)]=-z;
                        }
                    nn-=2;
                    break; // esce dalle its
                    }
                else {
                    if(its == 3*hqr_itsmax) {
                       alert("Mancata convergenza autovalore= "+nn+
                             "\niterazioni= "+its);
                        return -nn;
                        }
                    if(its== hqr_itsmax || its== 2*hqr_itsmax ) {
                         // shift eccezionale...
                         t=t+x;
                         for(i=1; i<=nn;i++) {
                             a.a[i+ld*i] -= x;
                             }
                         s= Math.abs(a.a[nn+ld*(nn-1)])+ 
                            Math.abs(a.a[nn-1+ld*(nn-2)]);
                         x=0.75*s;
                         y=x;
                         w=-0.4375*s*s;
                         }
                    for(j=nn-2; j>=l;j--) {
                        m=j;
                        z=a.a[m+ld*m];
                        r=x-z;
                        s=y-z;
                        p=(r*s-w)/a.a[m+1+ld*m]+a.a[m+ld*(m+1)];
                        q=a.a[m+1+ld*(m+1)]-z-r-s;
                        r=a.a[m+2+ld*(m+1)];
                        s=Math.abs(p)+Math.abs(q)+Math.abs(r);
                        p=p/s;
                        q=q/s;
                        r=r/s;
                        if(m==l) break;
                        u=Math.abs(a.a[m+ld*(m-1)])*(Math.abs(q)+Math.abs(r));
                        v=Math.abs(p)*(Math.abs(a.a[m-1+ld*(m-1)]) +
                          Math.abs(z)+Math.abs(a.a[m+1+ld*(m+1)]));
                        if(u+v<=v) break;
                        }
                    // 4
                    for(i=m+2;i<=nn;i++){
                        a.a[i+ld*(i-2)]=0.0;
                        if(i!=m+2) a.a[i+ld*(i-3)]=0.0;
                        }
                    for(k=m;k<nn;k++){
                        if(k!=m) {
                            p=a.a[k+ld*(k-1)];
                            q=a.a[k+1+ld*(k-1)];
                            r=0.0;
                            if(k!=nn-1) r=a.a[k+2+ld*(k-1)];
                            x=Math.abs(p)+Math.abs(q)+Math.abs(r);
                            if(x > 0.0) {
                                p=p/x;
                                q=q/x;
                                r=r/x;
                                }
                            }
                        s=Math.sqrt(p*p+q*q+r*r);
                        if(p<0.0) s=-s;
                        if(Math.abs(s)>0.0) {
                            if(k==m) {
                                if(l!=m) a.a[k+ld*(k-1)]=-a.a[k+ld*(k-1)];
                                }
                            else {
                                a.a[k+ld*(k-1)]=-s*x;
                                }
                            p=p+s;
                            x=p/s;
                            y=q/s;
                            z=r/s;
                            q=q/p;
                            r=r/p;
                            for(j=k;j<=nn;j++){ // Modifica le righe
                                p=a.a[k+ld*j]+q*a.a[k+1+ld*j];
                                if(k != nn-1) {
                                    p += r*a.a[k+2+ld*j];
                                    a.a[k+2+ld*j] -= p*z;
                                    }
                                a.a[k+1+ld*j] -= p*y;
                                a.a[k+ld*j] -= p*x;
                                }
                            j=k+3;
                            if(nn<j) j=nn;
                            for(i=l;i<=j;i++){  // Modifica le colonne
                                p=x*a.a[i+ld*k]+y*a.a[i+ld*(k+1)];
                                if(k != nn-1) {
                                    p += z*a.a[i+ld*(k+2)];
                                    a.a[i+ld*(k+2)] -= p*r;
                                    }
                                a.a[i+ld*(k+1)] -= p*q;
                                a.a[i+ld*k] -= p;
                                }
                            }
                        } // end for
                    } //end else
                } // end else
            } // end for its
        } // end for eigen
    } // end function
//
function iper(a){
    //
    // Scambia tra loro gli elementi dei due vettori di servizio
    // della matrice a. Serve per poter estrarre il vettore
    // della parte immaginaria degli autovalori complessi calcolati
    // dalla hqr.
    //
    var ael,j,ld=a.ld+1;
    for(j=1;j<=a.ndim;j++){
        ael=a.a[j];
        a.a[j]=a.a[j*ld];
        a.a[j*ld]=ael;
        }
    return a.ndim;
    }
//
var maxitera=50;          // usato in jacobi_eigen
var precisione=1.0e-14;   // usato in jacobi_eigen
//
function jacobi_eigen(a,v){
    //
    // Tradotto in JavaScript da Fortran pag. 460 di:
    // "Numerical Recipes in Fortran - The Art of Scientific Computing"
    // Second Edition (1992). 
    // W.H.Press, S.A.Teukolsky, W.T.Vetterling, B.P.Flannery.
    // ISBN 0-521-43064  Cambridge University Press.
    //
    // a e v debbono essere oggetti di tipo ARMA.
    // Il risultato e' il numero totale di rotazioni effettuate.
    // a.a in ingresso contiene la matrice simmetrica ed
    // in uscita contiene valori non significativi. 
    // a.a colonna 0, in uscita contiene gli eigenvalue NON ordinati.
    // v.a in uscita rappresenta la matrice degli autovettori.
    // Usa come ordine di a il valore di a.ndim.
    // Non supera il massimo di iterazioni maxitera.
    // ed usa uno zero efficace dato da precisione.
    //
    var i,j,n,lda,ldv,ip,iq,nrot,itera;
    var sm,tresh,g,h,t,theta,c,s,tau;
    n=a.ndim;
    ldv=v.ld+1;
    lda=a.ld+1;
    if((n+1)*(n+1) > v.a.length) {
        alert("Errore: la matrice destinata agli autovettori è troppo piccola \n"+
              "In jacobi_eigen : v.a.length ="+v.a.length+" Ordine di a ="+n);
        v.error=-1;
        return v.error;
        } 
    for(ip=1;ip<=n;ip++){
        for (iq=1;iq<=n;iq++){
            v.a[ip+iq*ldv]=0.0;
            }
        v.a[ip+ip*ldv]=1.0;
        }
    for(ip=1;ip<=n;ip++){
        a.a[ip*lda]=a.a[ip+ip*lda];
        a.a[ip]=a.a[ip*lda];
        v.a[ip]=0.0;
        }
    nrot=0;
    for (itera=0;itera<maxitera;itera++){
        // somma di tutti gli elementi sovradiagonali
        sm=0.0;
        for(ip=1; ip<n;ip++){
            for (iq=ip+1; iq <=n; iq++){
                sm=sm+Math.abs(a.a[ip+iq*lda])
                }
            }
        // alert("ip= "+ip+" iq= "+iq+" sm= "+sm);
        if(sm<precisione) {
        // alert("Somma dei sovradiagonali"+sm+"\nOrdine "+n+"; Iterazioni"+itera+"\n"+a.a);
            return nrot;
            } 
        tresh=0.0;
        if(itera < 4) tresh=0.2*sm/(n*n)
        //
        //  Per ogni elemento della triangolare superiore
        //
        for(ip=1;ip<n;ip++){
            for(iq=ip+1;iq<=n; iq++){
                g=100.0*Math.abs(a.a[ip+iq*lda]);
                if(itera>4 && 
                    (Math.abs(a.a[ip])+g <= Math.abs(a.a[ip])) &&
                    (Math.abs(a.a[iq])+g <= Math.abs(a.a[iq])) ) {
                    a.a[ip+iq*lda]=0.0;
                    continue;
                    }
                if(Math.abs(a.a[ip+iq*lda]) <= tresh) continue;
                // 
                //  A questo punto e' necessario effettuare una
                //  rotazione di Jacobi
                //
                h=a.a[iq]-a.a[ip];        
                if(Math.abs(h)+g<=Math.abs(h)) {
                    t=a.a[ip+iq*lda]/h;
                    }
                else {
                    theta=0.5*h/a.a[ip+iq*lda];
                    t=1.0/(Math.abs(theta)+Math.sqrt(1.0+theta*theta));
                    if(theta<0.0) t=-t;
                    }
                c=1.0/Math.sqrt(1.0+t*t);
                s=t*c;
                tau=s/(1.0+c);
                h=t*a.a[ip+iq*lda];
                v.a[ip]=v.a[ip]-h;
                v.a[iq]=v.a[iq]+h;
                a.a[ip]=a.a[ip]-h;
                a.a[iq]=a.a[iq]+h;
                a.a[ip+iq*lda]=0.0;
                for(j=1;j<ip;j++){
                    g=a.a[j+ip*lda];
                    h=a.a[j+iq*lda];
                    a.a[j+ip*lda]=g-s*(h+g*tau);
                    a.a[j+iq*lda]=h+s*(g-h*tau);
                    }
                for(j=ip+1;j<iq;j++){
                    g=a.a[ip+j*lda];
                    h=a.a[j+iq*lda];
                    a.a[ip+j*lda]=g-s*(h+g*tau);
                    a.a[j+iq*lda]=h+s*(g-h*tau);
                    }
                for(j=iq+1;j<=n;j++){
                    g=a.a[ip+j*lda];
                    h=a.a[iq+j*lda];
                    a.a[ip+j*lda]=g-s*(h+g*tau);
                    a.a[iq+j*lda]=h+s*(g-h*tau);
                    }
                for(j=1;j<=n;j++){
                    g=v.a[j+ip*ldv];
                    h=v.a[j+iq*ldv];
                    v.a[j+ip*ldv]=g-s*(h+g*tau);
                    v.a[j+iq*ldv]=h+s*(g-h*tau);
                    }
                nrot++;
                } // end iq
            } // end ip
            for(ip=1;ip<=n;ip++){
                a.a[ip*lda]=a.a[ip*lda]+v.a[ip];
                a.a[ip]=a.a[ip*lda];
                v.a[ip]=0.0;
                }
        } // end itera
        // Se arriva qui ha superato il numero limite di iterazioni
        alert("Raggiunto il limite di iterazioni "+maxitera);
    }
//
function lu_fa(a){
    //
    // NON usa i limiti di clip
    // Fattorizza a considerata di ordine a.ndim.
    //
    var j,k,n,ld,deter,t,ajj,dva,dtt;
    n=a.ndim;
    ld=a.ld+1;
    deter=1.0;
    for(j=1;j<=n;j++){
        for(jd=1; jd<j; jd++){
            for(id=j; id<=n; id++){
                a.a[id+j*ld]+=a.a[jd+j*ld]*a.a[id+jd*ld];
                }
            }
        dva=a.a[j+j*ld];
        dtt=dva*dva;
        k=j;
        for(i=j+1; i<=n; i++){
            dva=a.a[i+j*ld];
            dva=dva*dva;
            if(dva > dtt) {
                dtt=dva;
                k=i;
                }
            }
        a.a[j*ld]=k;
        if(0.0 >= dtt) {
            a.error=-j;
            a.scalare=0.0;
            alert("Esce per errore: dtt="+dtt);
            return -j;
            }
        if(k != j ){
            for (i=1; i<=n; i++){
                t=a.a[j+i*ld];
                a.a[j+i*ld]=a.a[k+i*ld];
                a.a[k+i*ld]=t;
                } 
            deter=-deter;
            }
        ajj=a.a[j+j*ld];
        deter=deter*ajj;
        if(0.0 >= Math.abs(deter)) {
            a.error=-j;
            a.scalare=deter;
            alert("Esce per errore: deter="+deter);
            return -j;
            }
        if(Math.abs(ajj) > 0.0) {
            ajj=1.0/ajj;
            a.a[j+j*ld]=ajj;
            }
        else {
            a.error=-j;
            a.scalare=deter;
            alert("Esce per errore: ajj="+ajj);
            return -j;
            }
        for (jd=1; jd<=j-1;jd++){
            for(id=j+1;id<=n;id++){
                a.a[j+id*ld]+=a.a[j+jd*ld]*a.a[jd+id*ld];
                }
            }
        ajj=-ajj;
        for (i=j+1; i<=n; i++){
            a.a[j+i*ld]=ajj*a.a[j+i*ld];
            } 
        } // end j
    a.scalare=deter;
    return 0;
    }
//
function lu_so(a){
    //
    // NON usa i limiti di clip
    // Il vettore termine noto viene sostituito dal vettore solunzione.
    //
    var j,k,n=a.ndim;
    var ld=a.ld+1;
    var xk;
    for (k=1; k <=n; k++){
        j=a.a[k*ld];
        xk=a.a[j];
        a.a[j]=a.a[k];
        a.a[k]=xk;
        }
    for (k=1; k <=n; k++){
        xk=-a.a[k]*a.a[k+k*ld];
        for(j=k+1; j <=n; j++){
            a.a[j]+=a.a[j+k*ld]*xk; 
            }
        a.a[k]=-xk;
        }
    for (k=n; k>=1; k--) {
        xk=a.a[k];
        for(j=1; j<=k-1;j++){
            a.a[j]+=a.a[j+k*ld]*xk;
            }
        }
    return 0;
    }
//
function matvet(a,v){
    //
    // NON usa i limiti di clip ma va da 1 al minore tra a.ndim e v.ndim e
    // sulle a.ld  righe di a. Prodotto matrice per vettore. 
    // Non azzera il vettore ma aggiunge.
    // Come vettore usa il vettore interno di v.
    // Come matrice usa la matrice a.
    // Il risultato viene aggiunto al vettore interno di a ossia
    // tra 1 e a.ld del vettore a.a
    //
    var j,k,i,h,m,vv;
    m=a.ndim;
    if(v.ndim < m) m=v.ndim;
    i=a.ld;
    for (j=1;j<=m;j++){
        i++;
        h=0;
        vv=v.a[j];
        for(k=1;k<=a.ld;k++){
            i++;
            h++;
            a.a[h]+=a.a[i]*vv;       
            }
        }
    a.error=0;
    }
//
function setallmat(a,s){
    //
    // NON usa i limiti di clip. Usa la clip_setmat su tutta la matrice.
    // Usa la dimensione di servizio a.ndim che non puo' superare
    // il valore di a.ordine.
    // Ricordarsi che s è una stringa contenente i dati.
    //
    var ss;
    clip(1,1,a.ndim,a.ndim)
    ss=clip_setmat(a,s);
    clip_old(a);
    return ss;
    }
//
function setallvet(v,s){
    //
    // NON usa i limiti di clip. Inizializza il vettore interno
    // usando la clip_setmat.
    // Ricordare che s è una stringa contenente i dati.
    //
    var ss;
    clip(v,1,0,v.ndim,0)
    ss=clip_setmat(v,s);
    clip_old(v)
    return ss;
    }
//
function setdim(a,n){
    //
    // Imposta o trova la dimensione di servizio usata quando
    // NON usa i limiti di clip. La dimensione di servizio non puo'
    // essere inferiore ad 1 o superiore ad a.ordine.
    //
    var k=parseInt(n)
    if(isNaN(k)) return a.ndim;
    a.ndim=k;
    if(a.ndim<=1) a.ndim=1;
    if(a.ndim>a.ordine) a.ndim=a.ordine;
    return k;
    }
//
function setvet(v,s,primo,ultimo){
    //
    // NON usa i limiti di clip. Inizializza il vettore interno da
    // primo a ultimo usando la clip_setmat.
    //
    var ss;
    clip(v,primo,0,ultimo,0)
    ss=clip_setmat(v,s);
    clip_old(a);
    return ss;
    }
//
function zerovet(v){
    //
    // NON usa i limiti di clip. Azzera l'intero vettore interno.
    //
    var j;
    for(j=1; j<=v.ndim;j++){
        v.a[j]=0.0;
        }
    }
//
// ____________________________________________________________
//
// Includo una piccola libreria pensata per operare
// e visualizzare vettori di vettori. 
// Il legame tra le due librerie è dato dalle funzioni
// davetdivet e invetdivet.
//
// LIBRERIA ELEMENTARE DI ALGEBRA LINEARE
//
// Questa libreria richiede che le matrici siano espresse come 
// vettori di vettori.
// Una versione piu' raffinata di libreria di algebra lineare
// sta nel file alin_libesterna(operativa).txt
//
function proscal(a,b){ // Prodotto scalare
   var j,n,ps;
   n=a.length;
   if(n>b.length) n=b.length;
   ps=0.0;
   for(j=0;n>j;j++){
       ps+=a[j]*b[j];
       }
   return ps;
   }
//
function matrixordine(a){  // Il minore delle lunghezze dei vettori 
   var j,n,m;
   m=a.length;
   n=m;
   for(j=0;m>j;j++){
       if(n>a[j].length) n=a[j].length;
       }
   return n;
   }
//
function matrixbo(c,x){ // matrice unitaria di Bottoni
   var i,j,jj,k,m,n,r,rn,a;
   m=x.length;
   n=matrixordine(c);
   if(m>n) m=n;
   k=0;
   for(j=0;m>j;j++){
       if(Math.abs(x[j])>Math.abs(x[k])) k=j;
       }
   r=x[k]*x[k];
   for(j=1;m>j;j++){
      jj=(j+k)%m;
      rn=r+x[jj]*x[jj];
      a=x[jj]/Math.sqrt(r*rn);
      for(i=0;m>i;i++)c[j][i]=0;
      for(i=0;j>i;i++){
          ii=(k+i)%m;
          c[j][ii]=x[ii]*a;
          }
      c[j][jj]=-Math.sqrt(r/rn); 
      r=rn;
      }
   a=1.0/Math.sqrt(r);
   for(i=0;m>i;i++)c[0][i]=a*x[i];
   }
//
function matrixho(c,x){ // matrice unitaria di Householder
   var m,n,r,d,dd,sg,xv,j,k;
   m=x.length;
   n=matrixordine(c);
   if(m>n) m=n;
   if(m==1){c[0][0]=1;return m;}
   r=0;
   for(j=0;m>j;j++){
       r+=x[j]*x[j]; }
   r=Math.sqrt(r);
   sg=1;
   if(0>x[0]) sg=-1;        
   d=sg/(r*(r+Math.abs(x[0])));
   for(j=1;m>j;j++){
       dd=d*x[j];
       for(k=1;j>=k;k++){
           c[j][k]=dd*x[k];
           c[k][j]=c[j][k];
           }
       c[j][j]-=sg;
       }
   xv=x[0]+sg*r;
   dd=xv*d;
   for(j=1;m>j;j++){
      c[0][j]=dd*x[j];
      c[j][0]=c[0][j]; }
   c[0][0]=xv*dd-sg;
   return m;
   }
//
function matrixprodotto(c,a,b){ // prodotto di matrici scritte come vettori di vettori
   var j,k,h,n,m,s;
   m=matrixordine(a);
   n=matrixordine(b);
   if(m>n) m=n;
   n=matrixordine(c);
   if(m>n) m=n;
   for(j=0;m>j;j++){
      for(k=0;m>k;k++){
          s=0.0;
          for(h=0;m>h;h++){
              s+=a[j][h]*b[h][k];
              }
          c[j][k]=s;
          }
      }
   return m;
   }
//
function matrixpervet(y,c,x){ // matrice per vettore 
   var m,n,j,k,s;
   m=x.length;
   n=y.length;
   if(m>n) m=n;
   n=matrixordine(c);
   if(m>n) m=n;
   for(j=0;m>j;j++){
       s=0;
       for(k=0;m>k;k++){
           s+=c[j][k]*x[k];
           }
       y[j]=s;
       }
   return m;
   }
//
function matrixtrasposta(c,a){ // trasposta
   var j,k,n,m;
   m=matrixordine(a);
   n=matrixordine(c);
   if(m>n) m=n;
   for(j=0;m>j;j++){
      for(k=0;m>k;k++){
          c[j][k]=a[k][j];
          }
      }
   return m;
   }
//
//  Per visualizzare i vettori di vettori
//
function tondo(a){ // arrotonda un numero
   var b,c,xtondo=5.0e-12;
   c=a;
   if(a>0){
       b=a+xtondo;
       if(1.e-11>b%1)c=b-b%1;
       }
   else {
       b=-a+xtondo;
       if(xtondo+xtondo>b%1)c=-(b-b%1);
       }
   return c;
   }
//
function matrixstampo(a){ // stampa una matrice quadrata
   var j,m,s;
   m=matrixordine(a);
   s="<table>"
   for(j=0;m>j;j++){
       if(j==0){
          s="<table><tr><td rowspan='"+m+"' align='center'>&#x250C;<br/>";
          for(k=0;m>k;k++){
              s+="&#x2502;<br/>"
              }
          s+="&#x2514;</td>";
          }
       else {
          s+="<tr>"
          } 
       for(k=0;m>k;k++){
           s+="<td align='center'>"+tondo(a[j][k])+"</td>";
           }
       if(j==0){
           s+="<td rowspan='"+m+"' align='center'>&#x2510;<br/>";
           for(k=0;m>k;k++){
               s+="&#x2502;<br/>"
               }
               s+="&#x2518;</td></tr>";
           }
       else {
           s+="</tr>";
           }
       }
   s+="</table>";
   return s;
   }
//
function vetstampo(v){ // stampa un vettore
   var m,s
   m=v.length;
   s="<br/>[ ";
   for(j=0;m>j;j++){
       s+=tondo(v[j])+", ";
       }
   s+="] ";
   return s;
   }
//
// ____________________________________________________________
//
// Alcune variabili comode per evitare ( a volte) di dover usare ![CDATA[
//
var stag=String.fromCharCode(60);
var etag=String.fromCharCode(62);
var stable=stag+"table border='1'"+etag;
var etable=stag+"/table"+etag;
var sdiv=stag+"div style='color:red'"+etag;
var ediv=stag+"/div"+etag;
//
// FINE LIBRERIA ALIN ( Algebra LINeare basata sull'oggetto ARMA ) ]]> 
//