//Content-type: application/x-javascript

//CONSTANTS
TOL=1e-14;

// MATH LIBRARIES //
N=function(x) {return parseFloat(x)}; //make strings numbers
nan2zero=function(x) {if (isNaN(x)) {return N(0)} else {return x}}; //make NaNs zeros
pol=function(x) {if(Math.abs(N(x))<TOL){return 0.0} else{return N(x)};}; //make smalls zeros

function norm(a) {
  //calculates the 2-norm of linear array a

  //first find the largest-magnitude element
  var index=-1;var si=0; var dum;
  for(i=0;i<a.length;i++) {
    if((dum=Math.abs(a[i]))>si) {
      si=dum; index=i;
    }
  }
  //si is now the magnitude of the largest-magnitude element
  if(si==0) return 0;
  else {
   var norm2 = 0;
   for(i=0;i<a.length;i++) {norm2 +=(a[i]/si)*(a[i]/si)};
   return si*Math.sqrt(norm2);
  }
}
function transpose(M) {
 //calculates the transpose of mxn matrix M
 var m=M.length;
 var n
 if(m>0) {n=M[0].length;}
 else    {n=0;}
 Mt=new Array(m);  
 for(i=0;i<m;i++) Mt[i]=new Array(n);
 for(var i=0;i<m;i++) {
  for(var j=0;j<n;j++) {
    Mt[j][i]=M[i][j];
  }
 }
 return Mt;
}
function trace(M) {
 //calculates the trace of mxm matrix M
 var m=M.length;
 var t=0;
 if(m>0) {
   for(var i=0;i<m;i++) {
     t+=M[i][i];
   }
 }
 return t;
}
function mxm(M1,M2) {
 //multiplies mxp matrix M1 with pxn matrix M2 to get mxn matrix M3
 var m=M1.length
 var p=M2.length
 if(p>0) {n=M2[0].length;}
 else    {n=0;}
 var M3=new Array(m); 
 for(i=0;i<m;i++) M3[i]=new Array(n);
 for(i=0;i<m;i++) {
   for(j=0;j<n;j++) {
     M3[i][j]=0.0;
     for(k=0;k<p;k++) {
       M3[i][j]+=M1[i][k]*M2[k][j];
     }
   }
 }
 return M3; 
}
function mpm(M1,M2) {
 //adds mxn matrix M1 with mxn matrix M2 to get mxn matrix M3
 var m=M1.length
 if(m>0) {n=M1[0].length;}
 else    {n=0;}
 var M3=new Array(m);
 for(i=0;i<m;i++) M3[i]=new Array(n);
 for(i=0;i<m;i++) {
   for(j=0;j<n;j++) {
     M3[i][j]=M1[i][j]+M2[i][j];
   }
 }
 return M3;
}
function sxm(s,M1) {
 //multiplies scalar s and mxn matrix M1 to get mxn matrix M3
 var m=M1.length
 if(m>0) {n=M1[0].length;}
 else    {n=0;}
 var M2=new Array(m);
 for(i=0;i<m;i++) M2[i]=new Array(n);
 for(i=0;i<m;i++) {
   for(j=0;j<n;j++) {
     M2[i][j]=s*M1[i][j];
   }
 }
 return M2;
}
function flatten(v) {
 //returns a linear array that corresponds to column vector v
 var l=v.length;
 a=new Array(l)
 for(i=0;i<l;i++) a[i]=v[i][0];
 return a;
}
function hat(v) {
 //returns a unit vector in the same direction as v, or zero if v=0
 var a=flatten(v);
 var n=norm(a);
 if(n>0) return sxm(1/n,v);
 else    return v;
}