function trace(v)
{
   document.write(v + "<br>");
}

function is_zero(v)
{
// return true if number is zero (i.e., close to zero)

   if (Math.abs(v) < 0.0000001)
      return true;
   else
      return false;
}

function is_equal(a,b)
{
// return true if both numbers are equal (within tolerance)

   return is_zero(a-b);
}

function is_zero_vector(v)
{
// return true if 2D vector is zero

   if (is_zero(v[0]) && is_zero(v[1]) )
      return true;

   return false;
}

function pseudoinv(v)
{
// return inverse or zero (when value is zero)

   if (is_zero(v))
      return 0.0;
   else
      return 1.0 / v;
}

function vector(x , y)
{
// create a vector
   var p = new Array(2);
   
   p[0] = x;
   p[1] = y;
   
   return p;
}

function vector_normalize(V)
{
// normalize 2D vector to unit length

   var m= Math.sqrt(V[0]*V[0] + V[1]*V[1]);

   if (!is_zero(m))
   {
      V[0] /= m;
      V[1] /= m;
   }
}

function vector_dot(v1, v2)
{
// vector dot product

   return v1[0]*v2[0] + v1[1]*v2[1];
}

function vector_magnitude(v)
{
// vector length

   return Math.sqrt(v[0]*v[0] + v[1]*v[1]);
}

function quadratic(a, b, c)
{
// solve quadratic equation (if a is zero, return [0,0])

   var X = new Array(2);

   X[0] = 0.0;
   X[1] = 0.0;

   if (!is_zero(a))
   { 
      var pm  = Math.sqrt(b*b - 4*a*c);
      var den = 1.0 / (2.0 * a);
      X[0] = (-b - pm) * den;
      X[1] = (-b + pm) * den;
   }

   return X;
}

function Matrix2x2(a00, a01, a10, a11)
{
// 2x2 matrix constructor

// allocate

   this.M    = new Array(2);
   this.M[0] = new Array(2);
   this.M[1] = new Array(2)

// initialize

   this.M[0][0] = a00;
   this.M[0][1] = a01;
   this.M[1][0] = a10;
   this.M[1][1] = a11;

   this.Set = function(a00, a01, a10, a11) 
   { 
      // setter method 

      this.M[0][0] = a00;
      this.M[0][1] = a01;
      this.M[1][0] = a10;
      this.M[1][1] = a11;
   }


   this.Determinant = function()
   {
      // return determinant of 2x2

      return this.M[0][0] * this.M[1][1] - this.M[0][1] * this.M[1][0];
   }

   this.Dump = function() 
   {
      // dump out the contents of this matrix

      document.write(this.M[0][0].toFixed(4) + " " + this.M[0][1].toFixed(4) + "<br>");
      document.write(this.M[1][0].toFixed(4) + " " + this.M[1][1].toFixed(4) + "<br><br>");
   }

   this.Transpose = function() 
   {
       // return transpose of this matrix

       return new Matrix2x2(this.M[0][0], this.M[1][0], this.M[0][1], this.M[1][1]);
   }

   this.Sub = function(B)
   {
      // return this matrix - B

      var C = new Matrix2x2(0, 0, 0, 0);

      C.M[0][0] = this.M[0][0] - B.M[0][0];
      C.M[0][1] = this.M[0][1] - B.M[0][1];
      C.M[1][0] = this.M[1][0] - B.M[1][0];
      C.M[1][1] = this.M[1][1] - B.M[1][1];
 
      return C;
   }

   this.Mul = function(B)
   {
      // return this matrix * B

      var C = new Matrix2x2(0, 0, 0, 0);

      C.M[0][0] = this.M[0][0] * B.M[0][0] + this.M[0][1]*B.M[1][0];
      C.M[0][1] = this.M[0][0] * B.M[0][1] + this.M[0][1]*B.M[1][1];
      C.M[1][0] = this.M[1][0] * B.M[0][0] + this.M[1][1]*B.M[1][0];
      C.M[1][1] = this.M[1][0] * B.M[0][1] + this.M[1][1]*B.M[1][1];
 
      return C;
   }

   this.IsEqual = function(B)
   {
      // return true if this matrix is equal to B

      if (!is_equal(this.M[0][0], B.M[0][0]))
         return false;

      return true;
   }

   this.IsIdentity = function()
   {
      // return true if this matrix is the identity matrix 

      return this.IsEqual(new Matrix2x2(1, 0, 0, 1));
   }

   this.IsOrthogonal = function()
   {
       // return true if this matrix is orthogonal

       return this.Mul(this.Transpose()).IsIdentity();
   }

   this.SwapCols = function()
   {
      // swap columns of this matrix

      var t = this.M[0][0];
      this.M[0][0] = this.M[0][1];
      this.M[0][1] = t;

      t = this.M[1][0];
      this.M[1][0] = this.M[1][1];
      this.M[1][1] = t;
   }

   this.SwapRows = function()
   {
      // swap rows of this matrix 

      var t = this.M[0][0];
      this.M[0][0] = this.M[1][0];
      this.M[1][0] = t;

      t = this.M[0][1];
      this.M[0][1] = this.M[1][1];
      this.M[1][1] = t;
   }

   this.Eigenvalues = function() 
   {
       // calculate eigenvalues of this matrix 
       // returns 2 element array of eigenvalues

       var za = 1.0;
       var zb = - (this.M[0][0] + this.M[1][1]);
       var zc = this.Determinant();

       return quadratic(za, zb, zc);
   }

   this.Eigenvector = function(lambda)
   {
      // given an eigenvalue, calculates corresponding 
      // eigenvector of this matrix

      // lambda * I - A = 0
      var E = new Matrix2x2(0,0,0,0);
      E.M[0][0] = lambda - this.M[0][0];
      E.M[0][1] = this.M[0][1];
      E.M[1][0] = this.M[1][0];
      E.M[1][1] = lambda - this.M[1][1];

      // normalize to unit vectors
      vector_normalize(E.M[0]);
      vector_normalize(E.M[1]);

      // allocate eigenvector
      var ev = new Array(2);

      if (!is_zero_vector(E.M[0]))
      {
          // return this non-zero vector
          ev[0] = E.M[0][1];
          ev[1] = E.M[0][0];
      }
      else if (!is_zero_vector(E.M[1]))
      { 
          // return this non-zero vector
          ev[0] = E.M[1][1];
          ev[1] = E.M[1][0];
      }
      else
      {
          // everything is zero
          ev[0] = 0.0;
          ev[1] = 0.0;
      }

      // return eigenvector
      return ev;
   }

   this.Transform = function(p)
   {
       // transform point
       
       // allocate q
       var q = new Array(2);
       
       // transform
       q[0] = this.M[0][0] * p[0] + this.M[0][1] * p[1];
       q[1] = this.M[1][0] * p[0] + this.M[1][1] * p[1];
       
       // return q
       return q;
   }

}

function SVD2x2(A)
{
   // calculate SVD of 2x2 matrix 'A'

   // get transpose of A
   var At  = A.Transpose();

   // calculate A*A'
   var AAt = A.Mul(At);

   // get eigenvalues of A*A'
   var eigenvalues = AAt.Eigenvalues();

   // this is the transpose of U (inverse too)
   // I'm calculating it first because the eigenvectors 
   // will be returned in the form of rows (which I can
   // simply assign as follows saving a little work).
   var Uinv = new Matrix2x2(0,0,0,0);
   Uinv.M[0] = AAt.Eigenvector(eigenvalues[0]);
   Uinv.M[1] = AAt.Eigenvector(eigenvalues[1]);

   // tranpose the eigenvectors in to the columns of U
   this.U = Uinv.Transpose();

   // calculate S (sigma) using the square roots of the eigenvalues of A*A'
   this.S = new Matrix2x2(Math.sqrt(eigenvalues[0]), 0, 0, Math.sqrt(eigenvalues[1]));

   // calculate pseudoinverse of S
   // if S is non-singular, Sinv is the inverse of S
   // if S is singular, Sinv is the pseudoinverse of S
   var Sinv = new Matrix2x2(pseudoinv(this.S.M[0][0]), 0, 0, pseudoinv(this.S.M[1][1]));

   // calculate V'
   // A = U*S*V'  -> inv(U)*A = inv(U)*U*S*V' -> inv(U)*A = S*V' ->
   // inv(S)*inv(U)*A = inv(S)*S*V' -> inv(S)*inv(U)*A = V'
   this.Vt = Sinv.Mul(Uinv).Mul(A);

   // SVD requires the diagonal of S in descending order
   // every swap on the diagonal requires a column swap in U
   // and a row swap in V'
   if (this.S.M[0][0] < this.S.M[1][1])
   {
      var t = this.S.M[0][0];
      this.S.M[0][0] = this.S.M[1][1];
      this.S.M[1][1] = t;
 
      this.U.SwapCols();
      this.Vt.SwapRows();
   }      

   // now let's try to clean up the ugly cases (when A is singular)
   if (is_zero_vector(this.Vt.M[0]))
   {
      // if first row of V' is zero, make it Orthogonal to the second row
      this.Vt.M[0][0] = -this.Vt.M[1][1];
      this.Vt.M[0][1] =  this.Vt.M[1][0];
   }

   if (is_zero_vector(this.Vt.M[1]))
   { 
      // if second row is zero, make it Orthogonal to the first row
      this.Vt.M[1][0] =  this.Vt.M[0][1];
      this.Vt.M[1][1] = -this.Vt.M[0][0];
   }

   if (is_zero_vector(this.Vt.M[0]) && is_zero_vector(this.Vt.M[1]))
   {
      // if either row of V is still zero (then they both are)
      // at this point, A is zero, so it doesn't matter
      // what V' is, so just set it to the standard basis / identity
      this.Vt.Set(1,0,0,1);
   }

   if (is_zero_vector(this.U.M[0]) && is_zero_vector(this.U.M[1]))
   {
      // if U is zero, A is zero, so U doesn't matter either,
      // so just set it to the standard basis / identity.
      this.U.Set(1, 0, 0, 1);
   }
}