/*
MatrixMath.js
C Odenthal
2008-04

This file contains JavaScript functions to do some matrix computations.
Nothing fancy, matrix multiplication and addition; maybe a bit more.

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at
your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT 
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License 
(at http://www.gnu.org/licences/lgpl.html) for more details.
*/

// Vectors are implemented as javascript Arrays of numbers. Vectors can
// be displayed as either row or column vectors. A matrix consists of an
// Array of row vectors. When a row or column is removed from a matrix
// it becomes nondenominational - either a row or a column vector, and
// a vector can be inserted into a matrix as either a row or column.
//
// All functions return explicit values and have no side effects. Since
// javascript passes Arrays by reference, explicit copies have to be
// made of all matrix and vector inputs to functions. This is inefficient,
// but we're not doing anything computationally intensive enough for this
// to matter. Much.

// Functions to create matrices and vectors.
// Many of these use the 'ran' function from ASCIIMathML.js

function rr () { return ran(-9,9,3) }

// The attempt is to use inheritance (rather than aggregation) from the
// builtin Array class to make Vector and Matrix classes. See page 116
// of "Javascript Objects" by Nakhimovsky and Myers.
//
// First some prototyped functions to add Array methods.
// They change the original Array and return the modified Array.
//

// Performs a random shuffle of the entries of the Array.
Array.prototype.shuffle = function () {
    if (this.copy) {
        that = this.copy()
    } else {
        that = this.slice()
    }
    for (var rnd, tmp, i=that.length-1; i>0; i--) {
        rnd = ran(0,i)
        tmp = that[rnd]
        that[rnd] = that[i]
        that[i] = tmp
    }
    return that
}

// Applies the provided function 'fun' to every entry of the Array.
// This prototype is provided by the Mozilla foundation and
// is distributed under the MIT license.
// http://www.ibiblio.org/pub/Linux/LICENSES/mit.license
// But I fiddled with the code's formatting. I couldn't help myself.
if (!Array.prototype.map) {
    Array.prototype.map = function(fun /*, thisp*/) {
        var len = this.length;
        if (typeof fun != "function")
            throw new TypeError();
        var res = new Array(len);
        var thisp = arguments[1];
        for (var i = 0; i < len; i++) {
            if (i in this)
                res[i] = fun.call(thisp, this[i], i, this);
        }
        return res;
    };
}

//
// Some numerical constants in my version of javascript:
// 1e-324   == 0
// 1e309    == Infinity
// 1+1e-15  == 1.000000000000001
// 1+1e-16  == 1
//

//
//  Vector constructor and methods
//

function Vector () {
    if (arguments.length==1) {
        for (var i=0; i<arguments[0].length; i++) {
            this.push(arguments[0][i])
        }
    }
}
Vector.prototype = new Array

//
// Methods that have side effects
// They change the original Vector and also return the modified Vector.
//

Vector.prototype.sortNum = function() {                     // Sorts the entries of a Vector
   return this.sort( function (a,b) { return a-b; } );      // numerically instead of using
}                                                           // the default alphabetical sorting.

Vector.prototype.extend = function (that) {
    for (var i=0; i<that.length; i++) {
        this.push(that[i])
    }
    return this
}

//
// Methods that don't have side effects
// They leave the original Vector alone and return a new Object.
//

Vector.prototype.copy = function () { return new Vector(this.slice()) }
Vector.prototype.diagonal = function (offset) {
    if (!offset) { offset = 0 }
    var m = this.length + Math.abs(offset)
    var mat = zeros(m,m)
    if (offset<0) {
        for (var i=0; i<this.length; i++) {
            mat[i-offset][i] = this[i]
        }
    } else {
        for (var i=0; i<this.length; i++) {
            mat[i][i+offset] = this[i]
        }
    }
    return mat
}

Vector.prototype.vandermonde = function (n) {
    if (!n) { n = this.length }
    var m = this.length
    var mat = constant_matrix(m,n,1)
    for (var i=0; i<m; i++) {
        for (var j=1; j<n; j++) {
            mat[i][j] = Math.pow(this[i],j)
        }
    }
    return mat
}

Vector.prototype.max = function () {
    var m = this[0]
    for (var i=0; i<this.length; i++) {
        m = Math.max(m,this[i])
    }
    return m
}

Vector.prototype.min = function () {
    var m = this[0]
    for (var i=0; i<this.length; i++) {
        m = Math.min(m,this[i])
    }
    return m
}

Vector.prototype.sum = function () {
    var s = 0
    for (var i=0; i<this.length; i++) {
        s += this[i]
    }
    return s
}

Vector.prototype.prod = function () {
    var p = 1
    for (var i=0; i<this.length; i++) {
        p *= this[i]
    }
    return p
}

Vector.prototype.add = function (that) {
    var vec = new Vector()
    if (!that.length) { 
        for (var i=0; i<this.length; i++) {
            vec.push(this[i]+that)
        }
    } else {
        for (var i=0; i<this.length; i++) {
            vec.push(this[i]+that[i])
        }
    }
    return vec
}

Vector.prototype.sub = function (that) {
    var vec = new Vector()
    if (!that.length) { 
        for (var i=0; i<this.length; i++) {
            vec.push(this[i]-that)
        }
    } else {
        for (var i=0; i<this.length; i++) {
            vec.push(this[i]-that[i])
        }
    }
    return vec
}

Vector.prototype.dot = function (that) {
    if (!that[0].length) { // row vector * col vector
        var sca = 0
        for (var i=0; i<this.length; i++) {
            sca += this[i]*that[i]
        }
        return sca
    } else { // row vector * matrix
        var vec = constant_vector(that[0].length,0)
        for (var i=0; i<this.length; i++) {
            vec = vec.add(that[i].mul(this[i]))
        }
        return vec
    }
}

Vector.prototype.mul = function (that) {
    if (!that.length) { // vector * scalar
        var vec = new Vector()
        for (var i=0; i<this.length; i++) {
            vec.push(that*this[i])
        }
    } else {
        if (!that[0].length) { // col vector * row vector
            mat = new Matrix()
            for (var i=0; i<this.length; i++) {
                mat.push(that.mul(this[i]))
            }
            return mat
        } else { // row vector * matrix
            var vec = constant_vector(that[0].length,0)
            for (var i=0; i<this.length; i++) {
                vec = vec.add(that[i].mul(this[i]))
            }
        }
    }
    return vec
}

Vector.prototype.div = function (that) {
    var vec = new Vector()
    if (!that.length) { 
        for (var i=0; i<this.length; i++) {
            vec.push(this[i]/that)
        }
    } else {
        for (var i=0; i<this.length; i++) {
            vec.push(this[i]/that[i])
        }
    }
    return vec
}

Vector.prototype.ascii = function (p) {
    return new Matrix([this]).transpose().ascii(p)
}

Vector.prototype.latex = function (p) {
    return new Matrix([this]).transpose().latex(p)
}

Vector.prototype.ascii_row = function (p) {
    return new Matrix([this]).ascii(p)
}

Vector.prototype.latex_row = function (p) {
    return new Matrix([this]).latex(p)
}

//
//  Matrix constructor and methods
//

function Matrix () {
    if (arguments.length==1) {
        for (var i=0; i<arguments[0].length; i++) {
            this.push(new Vector(arguments[0][i]))
        }
    }
}

Matrix.prototype = new Array

//
// Methods that have side effects. They change the original
// Matrix and also return it after modification.
//

//
// Methods that don't have side effects.
// They leave the original Matrix alone and return a new Object.
//

Matrix.prototype.stack = function (that) {
    // Leaves 'this' and 'that' unchanged.
    // Returns a matrix with the rows of 'this'
    // followed by the rows of 'that'.
    mat = this.copy()
    if (mat[0].length != that[0].length) {
        return new Matrix([[]])
    } else {
        that = new Matrix(that)
        for (var i=0; i<that.length; i++) {
            mat.push(that[i])
        }
        return mat
    }
}

Matrix.prototype.cat = function (that) {
    // Leaves 'this' and 'that' unchanged.
    // Returns a matrix with the columns of 'this'
    // followed by the columns of 'that'.
    mat = this.copy()
    if (mat.length != that.length) {
        return new Matrix([[]])
    } else {
        that = new Matrix(that)
        for (var i=0; i<mat.length; i++) {
            mat[i].extend(that[i])
        }
        return mat
    }
}

Matrix.prototype.tril = function (offset) {
    mat = this.copy()
    if (!offset) { offset = 0 }
    var m = mat.length
    var n = mat[0].length
    for (var i=0; i<m; i++) {
        for (var j=Math.max(0,1+offset+i); j<n; j++) {
            mat[i][j]=0
        }
    }
    return mat
}

Matrix.prototype.triu = function (offset) {
    mat = this.copy()
    if (!offset) { offset = 0 }
    var m = mat.length
    var n = mat[0].length
    for (var i=0; i<m; i++) {
        for (var j=0; j<Math.min(offset+i,n); j++) {
            mat[i][j]=0
        }
    }
    return mat
}

Matrix.prototype.augment = function (that) {
    var b = new Matrix([that]).transpose()
    return this.cat(b)
}

Matrix.prototype.flipud = function () {
    return this.copy().reverse()
}

Matrix.prototype.fliplr = function () {
    return this.transpose().reverse().transpose()
}

Matrix.prototype.mySlice = function (m1,m2,n1,n2) {
    mat = new Matrix()
    for (var i=m1; i<m2; i++) {
        mat.push(new Vector(this[i].slice(n1,n2)))
    }
    return mat
}

Matrix.prototype.copy = function () {
    return this.mySlice(0,this.length,0,this[0].length)
}

Matrix.prototype.diagonal = function (offset) {
    if (!offset) { offset = 0 }
    var vec = new Vector()
    if (offset<0) {
        n = Math.min(this.length + offset, this[0].length)
        for (var i=0; i<n; i++) {
            vec.push(this[i-offset][i])
        }
    } else {
        n = Math.min(this.length, this[0].length - offset)
        for (var i=0; i<n; i++) {
            vec.push(this[i][i+offset])
        }
    }
    return vec
}

Matrix.prototype.transpose = function () {
    var n = this.length
    var m = this[0].length
    var B = constant_matrix(m,n,0)
    for (var i=0; i<m; i++) {
        for (var j=0; j<n; j++) {
            B[i][j] = this[j][i]
        }
    }
    return B
}

Matrix.prototype.mul = function (that) {
    if (!that.length) { // matrix * scalar
        var mat = new Matrix()
        for (var i=0; i<this.length; i++) {
            mat.push(this[i].mul(that))
        }
        return mat
    } else { // matrix * col vector
        if (!that[0].length) {
            var vec = new Vector()
            for (var i=0; i<this.length; i++) {
                vec.push(this[i].dot(that))
            }
            return vec
        } else { // matrix * matrix
            var mat = new Matrix()
            for (var i=0; i<this.length; i++) {
                mat.push(this[i].mul(that))
            }
            return mat
        }
    }
}

Matrix.prototype.div = function (that) {
    if (!that.length) { // matrix / scalar
        var mat = new Matrix()
        for (var i=0; i<this.length; i++) {
            mat.push(this[i].div(that))
        }
        return mat
    } else { // matrix / col vector
        if (!that[0].length) {
            var vec = new Vector()
            for (var i=0; i<this.length; i++) {
                vec.push(this[i].div(that[i]))
            }
            return vec
        } else { // matrix * matrix
            var mat = new Matrix()
            for (var i=0; i<this.length; i++) {
                mat.push(this[i].div(that[i]))
            }
            return mat
        }
    }
}

Matrix.prototype.add = function (that) {
    var mat = new Matrix()
    if (!that.length) {
        for (var i=0; i<this.length; i++) {
            mat.push(this[i].add(that))
        }
    } else {
        for (var i=0; i<this.length; i++) {
            mat.push(this[i].add(that[i]))
        }
    }
    return mat
}

Matrix.prototype.sub = function (that) {
    var mat = new Matrix()
    if (!that.length) {
        for (var i=0; i<this.length; i++) {
            mat.push(this[i].sub(that))
        }
    } else {
        for (var i=0; i<this.length; i++) {
            mat.push(this[i].sub(that[i]))
        }
    }
    return mat
}

Matrix.prototype.isSquare = function () {
    return (this.length == this[0].length)
}

Matrix.prototype.isNonsingular = function () {
    if (this.isSquare()) {
        var R = this.mgs()[1]
        return R.isSquare()
    } else {
        return false
    }
}

Matrix.prototype.ldu = function () {
    var C = this.copy()
    var m = C.length
    var n = C[0].length
    var col = new Vector()
    var ddd = 0
    var row = new Vector()
    var Lt = new Matrix()
    var U = new Matrix()
    var D = new Vector()
    for (var i=0; i<n; i++) {
        // Here we look for a nonzero entry in column i.
        for (var k=0; k<m; k++) {
            // If we find one, use it and then go to the next column.
            if (C[k][i] != 0) {
                col = C.transpose()[i]
                row = C[k]
                ddd = row[i]
                Lt.push( col.div(ddd) )
                U.push(  row.div(ddd) )
                D.push( ddd )
                C = C.sub(col.mul(row.div(ddd)))
                break
            }
        }
    }
    return [Lt.transpose(), D.diagonal(), U]
}

Matrix.prototype.lu = function () {
    var C = this.copy()
    var m = C.length
    var n = C[0].length
    var col = new Vector()
    var ddd = 0
    var row = new Vector()
    var Lt = new Matrix()
    var U = new Matrix()

    for (var i=0; i<n; i++) {
        // Here we look for a nonzero entry in column i.
        for (var k=0; k<m; k++) {
            // If we find one, use it and then go to the next column.
            if (C[k][i] != 0) {
                col = C.transpose()[i]
                row = C[k]
                ddd = row[i]
                Lt.push( col.div(ddd) )
                U.push( row )

                C = C.sub(col.mul(row.div(ddd)))
                break
            }
        }
    }
    return [Lt.transpose(), U]
}

Matrix.prototype.gepp = function () {
    var C = this.copy()
    var m = C.length
    var n = C[0].length
    var col = new Vector()
    var ddd = 0
    var row = new Vector()
    var Lt = new Matrix()
    var U = new Matrix()
    var parity = 1
    var permuted = []

    for (var i=0; i<n; i++) {
        // Here we look for the largest entry in column i.
        var max_i = 0
        var max_v = 0
        for (var k=0; k<m; k++) {
            if (Math.abs(C[k][i]) > max_v) {
                max_i = k
                max_v = Math.abs(C[k][i])
            }
        }
        if (max_v > 0) {
            col = C.transpose()[i]
            row = C[max_i]
            ddd = row[i]
            col = col.div(ddd)

            Lt.push( col )
            U.push( row )
            C = C.sub(col.mul(row))

            if (i != max_i) {
                permuted.push(i)
                if (permuted.indexOf(max_i) == -1) {
                    parity = -parity
                    permuted.push(max_i)
                }
            }
        }
    }
    if (Lt.length == 0) {
        Lt = zeros(1,m)
        Lt[0][0] = 1
        U = zeros(1,n)
    }

    return [Lt.transpose(), U, parity]
}

Matrix.prototype.qdr = function () {
    // Returns [Q,D,R] where A=QDR and Q has orthogonal columns,
    // the diagonal D stores the inverses of the squares of the
    // column lengths of Q, and R is upper triangular.

    var C = this.copy()
    var m = C[0].length
    var Qt = new Matrix()
    var R = new Matrix()
    var row = new Vector()
    var col = new Vector()
    var d = new Vector()
    var sca = 0

    for (var i=0; i<m; i++) {
        // Here we look for a nonzero column.
        // If we find one, use it - but in any case go to the next column.
        col = C.transpose()[i]
        sca = col.dot(col)
        if (sca > 0.00000000000000000000000001) {
            row = col.mul(C)

            // Cleanup row.
            for (var j=0; j<i; j++) {
                row[j] = 0
            }
            // Store col, sca and row
            Qt.push( col )
            d.push( 1/sca )
            R.push( row )

            // Update C.
            C = C.sub(col.mul(row.div(sca)))
        }
    }
    return [Qt.transpose(), d.diagonal(), R]
}

Matrix.prototype.ratqdr = function () {
    // This assumes the original matrix has integer entries.
    // Returns [Q,D,R] where A=QDR and Q has orthogonal columns,
    // the diagonal D stores the inverses of the squares of the
    // column lengths of Q, and R is upper triangular.

    var C = this.copy()
    var m = C[0].length
    var Qt = new Matrix()
    var R = new Matrix()
    var row = new Vector()
    var col = new Vector()
    var d = new Vector()
    var sca = 0

    for (var i=0; i<m; i++) {
        // Here we look for a nonzero column.
        // If we find one, use it - but in any case go to the next column.
        col = C.transpose()[i]
        sca = col.dot(col)
        if (sca > 0.00000000000000000000000001) {
            row = col.mul(C)

            // Cleanup row.
            for (var j=0; j<i; j++) {
                row[j] = 0
            }
            // Store col, sca and row
            Qt.push( col )
            // d.push( 1/sca )
            d.push( new Fraction(1,sca) )
            R.push( row )

            // Update C.
            C = C.sub(col.mul(row.div(sca)))
        }
    }
    return [Qt.transpose(), d.diagonal(), R]
}

Matrix.prototype.qdir = function () {
    // Returns [Q,Di,R] where A=QDi^{-1}R and Q has orthogonal columns,
    // the diagonal Di stores the the squares of the
    // column lengths of Q, and R is upper triangular.

    var C = this.copy()
    var m = C[0].length
    var Qt = new Matrix()
    var R = new Matrix()
    var row = new Vector()
    var col = new Vector()
    var d = new Vector()
    var sca = 0

    for (var i=0; i<m; i++) {
        // Here we look for a nonzero column.
        // If we find one, use it - but in any case go to the next column.
        col = C.transpose()[i]
        sca = col.dot(col)
        if (sca > 0.00000000000000000000000001) {
            row = col.mul(C)

            // Cleanup row.
            for (var j=0; j<i; j++) {
                row[j] = 0
            }
            // Store col, sca and row
            Qt.push( col )
            d.push( sca )
            R.push( row )

            // Update C.
            C = C.sub(col.mul(row.div(sca)))
        }
    }
    return [Qt.transpose(), d.diagonal(), R]
}

Matrix.prototype.mgs = function () {

    var C = this.copy()
    var m = C[0].length
    var Qt = new Matrix()
    var R = new Matrix()
    var row = new Vector()
    var col = new Vector()
    var sca = 0

    for (var i=0; i<m; i++) {
        // Here we look for a nonzero column.
        // If we find one, use it - but in any case go to the next column.
        col = C.transpose()[i]
        sca = Math.sqrt( col.dot(col) )
        if (sca > 0.0000000000001) {
            row = col.mul(C)

            // scale row and column
            row = row.div(sca)
            col = col.div(sca)

            // Cleanup row.
            for (var j=0; j<i; j++) {
                row[j] = 0
            }
            // Store col and row.
            Qt.push( col )
            R.push( row )

            // Update C.
            C = C.sub(col.mul(row))
        }
    }
    return [Qt.transpose(), R]
}

Matrix.prototype.qr = function () {
    var Q = eye(this.length)
    var R = this.copy()
    var H = new Matrix()
    for (var j=0; j<R[0].length; j++) {
        H = R.householder(j)
        R = H.mul(R)
        Q = Q.mul(H)
    }
    return [Q,R]
}

Matrix.prototype.hessenberg = function () {
    var Q = eye(this.length)
    var T = this.copy()
    var H = new Matrix()
    for (var j=0; j<T.length-1; j++) {
        H = T.householder(j,j+1)
        T = H.mul(T.mul(H))
        Q = Q.mul(H)
    }
    return [Q,T]
}

Matrix.prototype.solve = function (b) { 

    if (b.length != this.length) { return false }

    if (!b[0].length) {
        var vec = true
        var B = new Matrix([b]).transpose()
    } else {
        var vec = false
        var B = new Matrix(b)
    }
    var N = B[0].length

    var RY = this.cat(B).mgs()[1]
    var R = RY.mySlice(0,RY.length,0,this[0].length)
    var r = R.diagonal()

    // return r

    if (r.prod() == 0) { return false }

    R = R.mySlice(0,r.length,0,r.length)
    var Y = RY.mySlice(0,r.length,R[0].length,RY[0].length)
    var X = R.backSub(Y)

    if (vec == true) { X = X.transpose()[0] }

    return X
}

Matrix.prototype.luSolve = function (b) { 

    if (!this.isSquare()) {
        var At = this.transpose()
        var AtA = At.mul(this)
        var Atb = At.mul(b)
        return AtA.luSolve(Atb)
    }

    var LU = this.lu()
    var L = LU[0]
    var U = LU[1]

    if (!U.isSquare()) { return this.qrSolve(b) }

    return U.backSub(L.forwardSub(b))
}

Matrix.prototype.qrSolve = function (b) { 

    var QR = this.mgs()
    var Q = QR[0]
    var R = QR[1]

    if (!R.isSquare()) { return false }

    var Qt = Q.transpose()
    var c = Qt.mul(b)

    return R.backSub(c)
}

Matrix.prototype.qdrSolve = function (b) { 

    var QDR = this.qdr()
    var Q = QDR[0]
    var R = QDR[2]

    if (!R.isSquare()) { return false }

    var Qt = Q.transpose()
    var c = Qt.mul(b)

    return R.backSub(c)
}

Matrix.prototype.backSub = function (b) {
    if (!b[0].length) {
        var x = new Vector(b)
    } else {
        var x = new Matrix(b)
    }

    // Back substitute while assuming that
    // 'this' is upper triangular and invertible.

    if (!x[0].length) {
        for (var k=this.length-1; k>-1; k--) {
            x[k] = x[k]/this[k][k]
            for (var i=0; i<k; i++) {
                x[i] = x[i] - x[k]*this[i][k]
            }
        }
    } else {
        for (var k=this.length-1; k>-1; k--) {
            x[k] = x[k].div(this[k][k])
            for (var i=0; i<k; i++) {
                x[i] = x[i].sub(x[k].mul(this[i][k]))
            }
        }
    }
    return x
}

Matrix.prototype.forwardSub = function (b) {
    // Forward substitute while assuming that
    // 'this' is lower triangular and invertible.

    // Make a Matrix U that has the rows and
    // columns of A reversed. So U=JAJ where
    // J is the identity matrix with rows reversed.
    // Then Ux=Jb means A(Jx) = b
    var U = new Matrix()
    for (var i=0; i<this.length; i++) {
        U.push(this[i].copy().reverse())
    }
    U.reverse()

    return U.backSub(b.copy().reverse()).reverse()
}

Matrix.prototype.inv = function () {
    if (!this.isSquare()) {
        return false
    } else {
        return this.luSolve(eye(this.length))
    }
}

Matrix.prototype.det = function () {
    if (!this.isSquare()) { return false }
    if (!this.isNonsingular()) { return 0 }

    // If we've gotten this far, 'this' is a square matrix.
    var LUp = this.gepp()

    var dL = LUp[2]          // determinant of L (1 or -1)

    var U   = LUp[1]        // U upper triangular
    var dU = U.diagonal().prod()

    if (!U.isSquare()) {
        return 0
    } else { 
        return dL*dU
    }
}

Matrix.prototype.intDet = function () {
    return Math.round(this.det())
}

Matrix.prototype.jacobi = function (i,j) {
    // If Matrix is a symmetric matrix, then
    // this returns a Jacobi rotation matrix J
    // such that the (i,j) and (j,i) terms of
    // J^T*'this'*J would be 0.

    var mat = eye(this.length)

    if (i != j) {
        var a = this[i][i]
        var d = this[i][j]
        var b = this[j][j]
        if (Math.abs(d) > 1e-25) {
            var tau = (a-b)/(2*d)
            var t = sign(tau)/( Math.abs(tau) + Math.sqrt(1+tau*tau) )
            var c = 1/Math.sqrt(1+t*t)
            var s = c*t

            mat[i][i] = c
            mat[i][j] = -s
            mat[j][i] = s
            mat[j][j] = c
        }
    }

    return mat
}

Matrix.prototype.householder = function (j,k) {
    // Returns a Householder reflector H such that
    // H*'this' has zeros in column j past position
    // k.

    if (!k) {
        k = j
    }

    var mat = eye(this.length)
    var vec = this.transpose()[j]
    for (var i=0; i<k; i++) {
        vec[i] = 0
    }
    var v = Math.sqrt(vec.dot(vec))
    if (v==0) { return mat }

    if (vec[k]<0) { v = -v }
    vec[k] += v

    // I-2(v v^t)/(v^t v)
    return mat.sub(vec.mul(vec.mul(2).div(vec.dot(vec))))

}

// Formatting methods for the Matrix class

// The next method takes a precision 'p' and produces a
// string that will be formatted by latex as a matrix with
// round brackets and (at most) 'p' digit fixed point
// numbers.

Matrix.prototype.platex = function (p) {
    if (!p) { p = 4 }
    var m = this.length
    var n = this[0].length
    var As = []
    var t = 0
    for (var i=0; i<m; i++) {
        var row = []
        for (var j=0; j<n; j++) {
            row.push( my_fixed(this[i][j],p) )
        }
        As.push(row.join( " & " ))
    }
    As = "\\begin{pmatrix} " + As.join( " \\\\ " ) + " \\end{pmatrix}"
    return As
}

Matrix.prototype.latex = function (aug,p) {
    if (!p) { p = 4 }
    var m = this.length
    var n = this[0].length
    var As = []
    var t = 0
    for (var i=0; i<m; i++) {
        var row = []
        for (var j=0; j<n; j++) {
            row.push( my_fixed(this[i][j],p) )
        }
        As.push(row.join( " & " ))
    }
    var align = ""
    for (var i=0; i<n-1; i++) {
        align += "r"
    }
    if (!aug) {
        align += "r"
    } else {
        align += "|r"
    }
    As = "\\left(\\begin{array}{" + align + "} " + As.join( " \\\\ " ) + " \\end{array}\\right)"
    return As
}

Matrix.prototype.natlatex = function (aug,p) {
    // if (!p) { p = 4 }
    var m = this.length
    var n = this[0].length
    var As = []
    var t = 0
    for (var i=0; i<m; i++) {
        var row = []
        for (var j=0; j<n; j++) {
            row.push( this[i][j] )
        }
        As.push(row.join( " & " ))
    }
    var align = ""
    for (var i=0; i<n-1; i++) {
        align += "r"
    }
    if (!aug) {
        align += "r"
    } else {
        align += "|r"
    }
    As = "\\left(\\begin{array}{" + align + "} " + As.join( " \\\\ " ) + " \\end{array}\\right)"
    return As
}

// The next method takes a precision 'p' and produces a
// string that will be formatted by ASCIIjs as a matrix with
// round brackets and (at most) 'p' digit fixed point
// numbers.

Matrix.prototype.ascii = function (p) {
    if (!p) { p = 4 }
    var m = this.length
    var n = this[0].length
    var As = []
    var t = 0
    for (var i=0; i<m; i++) {
        var row = []
        for (var j=0; j<n; j++) {
            row.push( my_fixed(this[i][j],p) )
        }
        As.push(row)
    }
    As = "(" + As.join( "),(" ) + ")"
    if (m>1) { As = "(" + As + ")" }
    return As
}

Matrix.prototype.natascii = function () {
    var m = this.length
    var n = this[0].length
    var As = []
    for (var i=0; i<m; i++) {
        var row = []
        for (var j=0; j<n; j++) {
            row.push( this[i][j] )
        }
        As.push(row)
    }
    As = "(" + As.join( "),(" ) + ")"
    if (m>1) { As = "(" + As + ")" }
    return As
}

// The next method produces a string representation that
// is a javascript Array of Arrays. ASCIIMathML.js will
// format this as a matrix with square brackets and
// floating point numbers.

Matrix.prototype.jscript = function () {
    var m = this.length
    var n = this[0].length
    var As = []
    var t = 0
    for (var i=0; i<m; i++) {
        var row = []
        for (var j=0; j<n; j++) {
            row.push( String(this[i][j]) )
        }
        As.push(row)
    }
    As = "[" + As.join( "],[" ) + "]"
    if (m>1) { As = "[" + As + "]" }
    return As
}

//
// Global functions using Matrix and Vector classes
//

function random_vector (m,lo,hi) {
    var vec = new Vector()
    for (var i=0; i<m; i++) {
        vec.push( ran(lo,hi) )
    }
    return vec
}

function random_matrix (m,n,lo,hi) {
    var mat = new Matrix()
    for (var i=0; i<m; i++) {
        mat.push( random_vector(n,lo,hi) )
    }
    return mat
}

function constant_vector (n,s) {
    var vec = new Vector()
    for (var i=0; i<n; i++) {
        vec.push(s)
    }
    return vec
}

function constant_hector (n,s,fraction) {
    var vec = new Hector([],fraction)
    s = new vec.f(s)
    for (var i=0; i<n; i++) {
        vec.push(s)
    }
    return vec
}

function constant_matrix (m,n,s) {
    var mat = new Matrix()
    for (var i=0; i<m ; i++) {
        mat.push( constant_vector(n,s) )
        }
    return mat
}

function pm_vector (n) {
    // vector with +1 and -1 entries
    var vec = new Vector()
    for (var i=0; i<n; i++) {
        vec.push(2*ran(0,1)-1)
    }
    return vec
}

function pm_matrix (m,n) {
    // matrix with +1 and -1 entries
    var mat = new Matrix()
    for (var i=0; i<m; i++) {
        mat.push(pm_vector(n))
    }
    return mat
}

function scalar_matrix (m,n,s) {
    var mat = constant_matrix (m,n,0)
    var k = Math.min(m,n)
    for (var i=0; i<k; i++) { mat[i][i] = s }
    return mat
}

function zeros (m,n) { return constant_matrix (m,n,0) }
function eye (m) { return scalar_matrix (m,m,1) }

function stochastic_matrix (n,hi) {
    if (!hi) { hi = 10 }
    var mat = new Matrix()
    function t (m) { return ran(1,Math.max(1,Math.abs(m))) }
    for (var i=0; i<n; i++) {
        var row = new Vector()
        var sum = 0
        for (var j=0; j<n; j++) {
            row.push(t(hi-sum))
            sum += row[j]
        }
        if (sum <= 10) {
            row[i] += (10 - sum)
            sum = 10
        } else {
            row[i] += (100 - sum)
            sum = 100
        }
        for (var j=0; j<n; j++) {
            row[j] /= sum
        }
        mat.push(row)
    }
    return mat.transpose()
}

function pivot_matrix (m,n) {
    var A = zeros (1,n)
    var j = 0
    var p = ran(0,2)
    while ((p == 0) && (j < n-1)) {
        p = ran(0,1)
        j += 1
    }
    if (p != 0) { A[0][j] = 1 }

    if (m == 1) { return A }

    // We reach here if m>1
    if (j == n-1) { return A.stack(zeros(m-1,n)) }

    // We reach here if m>1 and j<n-1
    var B = zeros(m-1,j+1).cat(pivot_matrix(m-1,n-j-1))
    return A.stack(B)
}

// This produces a lower triangular matrix with ones down the diagonal.
function lower_triangular_matrix (m,n,lo,hi) {
    var mat = random_matrix (m,n,lo,hi).tril()
    var r = Math.min(m,n)
    for (var i=0; i<r; i++) {
        mat[i][i] = 1
    }
    return mat
}

// This produces an upper triangular matrix with
// nonzero entries down the diagonal.
function upper_triangular_matrix (m,n,lo,hi) {
    var mat = random_matrix (m,n,lo,hi).triu()
    var r = Math.min(m,n)
    for (var i=0; i<r; i++) {
        while (mat[i][i] == 0) {
            mat[i][i] = ran(lo,hi)
        }
    }
    return mat
}

function known_rank_matrix (m,n,r,lo,hi) {
    var lmat = lower_triangular_matrix (m,r,lo,hi)
    var rmat = upper_triangular_matrix (r,n,lo,hi)
    return lmat.mul(rmat)
}

function nonsingular_matrix (m,lo,hi) {
    var mat = random_matrix (m,m,lo,hi)
    while (!mat.isNonsingular()) {
        mat = random_matrix (m,m,lo,hi)
    }
    return mat
}

function symmetric_matrix (m,lo,hi) {
    var mat = lower_triangular_matrix(m,m,lo,hi)
    mat = mat.add(mat.transpose())
    for (var i=0; i<m; i++) {
        mat[i][i] = ran(lo,hi)
    }
    return mat
}

function positive_definite_matrix (m,hi) {
    if (!hi) { hi = 2 }
    var M = nonsingular_matrix(m,m,-hi,hi)
    var D = random_vector(m,1,hi+1).diagonal()
    return M.mul(D.mul(M.transpose()))
}

function orthogonal_cols_matrix(m,n) {
    switch(n) {
        case 1:
            x = 9
            break

        case 2:
            x = 3
            break

        case 3:
            x = 2
            break

        default:
            x = 1
            break
    }
    var At = known_rank_matrix(n,m,n,-x,x)
    for (var i=0; i<n; i++) {
        for (var j=i+1; j<n; j++) {
            var a = At[i].dot(At[j])
            var b = At[i].dot(At[i])
            var d = gcd(a,b)
            At[j] = At[j].mul(b/d).sub(At[i].mul(a/d))
        }
    }

    for (var i=0; i<n; i++) {
        d = Agcd(At[i])
        At[i] = At[i].div(d)
    }

    return At.shuffle().transpose()
}

function ref_matrix (m,n,lo,hi) {
    var p_mat = pivot_matrix(m,n)

    var U = lower_triangular_matrix(n,n,lo,hi).transpose()
    var D = pm_vector(m).diagonal().mul(random_vector(m,1,hi)).diagonal()
    var A = D.mul(p_mat.mul(U))

    return A
}

function rref_matrix (m,n,lo,hi) {
    var p_mat = pivot_matrix(m,n)

    // We get the list of pivot columns.
    var p_col = []
    for (var i=0; i<m; i++) {
        for (var j=i; j<n; j++) {
            if (p_mat[i][j] == 1) {
                p_col.push(j)
                break
            }
        }
    }

    var U = lower_triangular_matrix(n,n,lo,hi).transpose()
    var A = p_mat.mul(U)

    // Put the zeros back in the pivot columns.
    for (var i=0; i<p_col.length; i++) {
        for (var k=0; k<i; k++) {
            A[k][p_col[i]] = 0
        }
    }

    return [A,p_col]
}

// This probably could use some work. Maybe later.
//
function my_calculate(inputId,outputId) {
  var str = document.getElementById(inputId).value;
  var err = "";
  var ind = str.lastIndexOf("\n");
  if (ind==str.length-1) str = str.slice(0,ind);
  str = str.slice(str.lastIndexOf("\n")+1);
  try {
    var res = eval(mathjs(str));
  } catch(e) {
    err = "syntax incomplete";
  }
  if (!isNaN(res) && res!="Infinity") 
    str = "`"+str+" ` "+(Math.abs(res-Math.round(res*1000000)/1000000)<1e-15?Math.round(res*1000000)/1000000:res)+err; 
  else if (str!="") str = "`"+str+"` = undefined"; //debug:+mathjs(str);
  var outnode = document.getElementById(outputId);
  var n = outnode.childNodes.length;
  for (var i=0; i<n; i++)
    outnode.removeChild(outnode.firstChild);
  outnode.appendChild(document.createTextNode(str));
  AMprocessNode(outnode);
}

//  The next function differs from the 'chop' function in ASCIIMathML.js
//  in that it rounds instead of truncating the number. It differs from
//  the built-in 'toFixed' method in that it drops trailing zeros.

function my_fixed (x,p) {
    var w = 1
    for (var i=0; i<p; i++) { w *= 10 }
    return Math.round(x*w) / w
}

function gcd (a,b) {
    if (b==0) {
        return Math.abs(a)
    } else {
        return gcd(b,a%b)
    }
}

function Agcd (A) {
    var d = A[0]
    for (i=1; i<A.length; i++) {
        d = gcd(d,A[i])
    }
    return d
}

function lcm (a,b) {
    a = Math.abs(a)
    b = Math.abs(b)
    return a*b/gcd(a,b)
}

///////////////////////////////////////////////////
// These two functions are redundant (I think) due
// to the 'sortNum' method for Vectors.

function sort_by_number_asc(a,b) {
    return a-b
}

function sort_by_number_desc(a,b) {
    return b-a
}


///////////////////////////////////////////////////
// Code for rational numbers
//

// The next function is redundant (I think) because
// it is also in the ASCIIMathML.js file.

var sign = function(x) { return (x==0?0:(x<0?-1:1)) };

function Fraction (numer,denom) {

    // make sure there is a denominator
    if (arguments.length == 1) {
        if (!numer.denominator) {
            denom = 1
        } else {
            denom = numer.denominator
            numer = numer.numerator
        }
    }

    // make denominator positive
    var s = sign(denom)
    numer = numer * s
    denom = denom * s

    // reduce fraction
    var g = gcd(numer,denom)
    numer = numer / g
    denom = denom / g

    this.numerator = numer
    this.denominator = denom
}

Fraction.prototype = new Object

Fraction.prototype.toString = function () {
    fraction_str = this.numerator.toString()
	if (this.denominator != 1) {
        fraction_str = fraction_str + "/" + this.denominator.toString()
    }
    return fraction_str
}

Fraction.prototype.valueOf = function () {
    return this.numerator / this.denominator
}

Fraction.prototype.add = function (other) {
    if (!other.denominator) {
        if (Math.round(other) != other) {
            return this.valueOf() + other
        } else {
            other = new Fraction(other)
        }
    }
    numer = this.numerator * other.denominator + this.denominator * other.numerator
    denom = this.denominator * other.denominator
    return new Fraction(numer,denom)
}

Fraction.prototype.neg = function () {
    numer = -this.numerator
    denom = this.denominator
    return new Fraction(numer,denom)
}

Fraction.prototype.sub = function (other) {
    if (!other.denominator) {
        if (Math.round(other) != other) {
            return this.valueOf() - other
        } else {
            other = new Fraction(other)
        }
    }
    numer = this.numerator * other.denominator - this.denominator * other.numerator
    denom = this.denominator * other.denominator
    return new Fraction(numer,denom)
}

Fraction.prototype.mul = function (other) {
    if (!other.denominator) {
        if (Math.round(other) != other) {
            return this.valueOf() * other
        } else {
            other = new Fraction(other)
        }
    }
    numer = this.numerator * other.numerator
    denom = this.denominator * other.denominator
    return new Fraction(numer,denom)
}

Fraction.prototype.inv = function () {
    numer = this.denominator
    denom = this.numerator
    return new Fraction(numer,denom)
}

Fraction.prototype.div = function (other) {
    if (!other.denominator) {
        if (Math.round(other) != other) {
            return this.valueOf() / other
        } else {
            other = new Fraction(other)
        }
    }
    numer = this.numerator * other.denominator
    denom = this.denominator * other.numerator
    return new Fraction(numer,denom)
}

Fraction.prototype.pow = function (expon) {
    var numer = Math.pow(this.numerator,expon)
    var denom = Math.pow(this.denominator,expon)
    return new Fraction(numer,denom)
}

///////////////////////////////////////////////////
// Code for decimal numbers
//

function Decimal (value,prec) {

    // make sure there is a precision
    if (arguments.length == 1) { prec = 4 }

    this.value = value
    this.prec = prec
}

Decimal.prototype = new Object

Decimal.prototype.toString = function () {
    var w = Math.pow(10,this.prec)
    var r = Math.round(this.value * w) / w
    return r.toString()
}

Decimal.prototype.valueOf = function () {
    return this.value
}

Decimal.prototype.add = function (other) {
    return new Decimal(this.value + other, this.prec)
}

Decimal.prototype.neg = function () {
    return new Decimal(-this.value, this.prec)
}

Decimal.prototype.sub = function (other) {
    return new Decimal(this.value - other, this.prec)
}

Decimal.prototype.mul = function (other) {
    return new Decimal(this.value * other, this.prec)
}

Decimal.prototype.inv = function () {
    return new Decimal(1 / this.value, this.prec)
}

Decimal.prototype.div = function (other) {
    return new Decimal(this.value / other, this.prec)
}

Decimal.prototype.pow = function (expon) {
    number = Math.pow(this, expon)
    return new Decimal(number, this.prec)
}