package jamaica;

import java.io.PrintWriter;
import java.io.RandomAccessFile;
import java.io.IOException;
import java.text.NumberFormat;
import java.text.DecimalFormat;
import java.text.DecimalFormatSymbols;
import java.util.Locale;

/**
 * Jamaica - (yet another) Java Matrix Class
 *
 * @author Hanhua Feng - hanhua@cs.columbia.edu
 * @version $Id: Matrix.java,v 1.27 2003/05/21 15:32:24 hanhua Exp $
 */

public class Matrix implements Cloneable {

    /* internal data structure
     * a:  array to store matrix entries
     * m:  number of rows
     * n:  number of columns
     * r:  index (in the array) of the first entry of the matrix
     * ms: index difference if the row number increments
     * ns: index difference if the column number increments
     */
    double [] a;
    int m, n, r, ms, ns;

    /** Construct an m-by-n matrix and fill with zeros.
     * @param m     number of rows
     * @param n     number of columns
     */
    public Matrix( int m, int n ) {
        this.m = m;
        this.n = n;
        r = 0;
        ms = 1;     // the matrix is saved column-by-column
        ns = m;
        a = new double[m*n];
    }

    /** Construct an m-by-n matrix and fill with a scalar value.
     * @param m     number of rows
     * @param n     number of columns
     * @param s     the scalar value to be filled
     */
    public Matrix( int m, int n, double s ) {
        this.m = m;
        this.n = n;
        r = 0;
        ms = 1;
        ns = m;
        a = new double[m*n];
        for ( int i=0; i<a.length; i++ )
            a[i] = s;
    }

    /** Construct an m-by-n matrix with an array .
     * @param m     number of rows
     * @param n     number of columns
     * @param array the initial array (contructor will not make copies).
     * @param byrow the array is stored by row (default: by column)
     */
    public Matrix( int m, int n, double[] array, boolean byrow ) {
        this.m = m;
        this.n = n;
        r = 0;
        if ( byrow )
        {
            ms = n;
            ns = 1;
        }
        else
        {
            ms = 1;
            ns = m;
        }
        a = array;
    }

    /** Construct a matrix sharing to the same data of another matrix.
     * @param x     original matrix
     */
    public Matrix( Matrix x ) {
        m = x.m;
        n = x.n;
        r = x.r;
        ms = x.ms;
        ns = x.ns;
        a = x.a;
    }

    /** Construct a matrix to be a submatrix of another matrix
     * @param x     original matrix
     * @param mr    row range 
     * @param nr    column range
     */
    public Matrix( Matrix x, Range mr, Range nr ) {
        if ( ! mr.in( 0, x.m-1 ) || ! nr.in( 0, x.n-1 ) )
        {
            throw new IllegalArgumentException( 
                "Invalid range: " + x.m + "*" + x.n 
                + " [" + mr.toString() + "," + nr.toString() + "]" );
        }

        a = x.a;
        m = mr.size();
        n = nr.size();
        r = x.r + mr.first() * x.ms + nr.first() * x.ns;
        ms = x.ms * mr.interval();
        ns = x.ns * nr.interval();
    }

    /** Return an n-by-n identity matrix
     * @param n     number of rows or columns
     * @return      n-by-n identity matrix
     */
    public static Matrix eye( int n ) {
        Matrix x = new Matrix( n, n );
        for ( int i=0; i<n; i++ )
            x.set( i, i, 1.0 );
        return x;
    }

    /** Return an m-by-n matrix with random entries
     * @param m     number of rows
     * @param n     nubmer of columns
     * @return      m-by-n random matrix
     */
    public static Matrix random( int m, int n ) {
        Matrix x = new Matrix( m, n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                x.set( i, j, Math.random() );
        return x;
    }

    /** Return an m-by-n matrix with incremental entries.  All
     * differences between pairs of two vertically adjacent entries
     * are the same, and so for for horizontally pairs.
     * @param m    number of rows
     * @param n    number of columns
     * @param s    the value of entry (0,0)
     * @param dm   vertical increment 
     * @param dn   horizontal increment
     * @return     new m-by-n matrix
     */
    public static Matrix indgen( int m, int n, double s, 
                                 double dm, double dn ) {
        Matrix x = new Matrix( m, n );
        for ( int j=0; j<n; j++ )
        {
            double p = s;
            for ( int i=0; i<m; i++, p+=dm )
                x.set( i, j, p );
            s += dn;
        }
        return x;
    }

    // Return an string indicating an error for array dimensions
    private final String dimErrMsg() {
        return "Invalid dimensions: " + m + "*" + n;
    }

    private final String dimErrMsg( Matrix b ) {
        return "Dimensions not match: " + m + "*" + n + " <=> " 
            + b.m + "*" + b.n;
    }

    // get correspoding bit in bmk for matrix element (i,j)
    final boolean getBitMask( BitArray bmk, int i, int j ) {
        return bmk.get( i+j*m );
    }

    // set correspoding bit in bmk for matrix element (i,j)
    private final void setBitMask( BitArray bmk, int i, int j ) {
        bmk.set( i+j*m );
    }

    /** Get row dimensions 
     * @return      the number of rows
     */
    public final int height() {
        return m;
    }

    /** Get column dimensions 
     * @return      the number of rows
     */
    public final int width() {
        return n;
    }

    /** Get an entry of the matrix
     * @param i     row index (zero-based)
     * @param j     column index (zero-based)
     */
    public final double get( int i, int j ) {
        return a[r+i*ms+j*ns];
    }

    /** Assign a value to an entry of the matrix
     * @param i     row index (zero-based)
     * @param j     column index (zero-based)
     * @param x     value to be assigned
     */
    public final void set( int i, int j, double x ) {
        a[r+i*ms+j*ns] = x;
    }

    /** Add a value to an entry of the matrix
     * @param i     row index (zero-based)
     * @param j     column index (zero-based)
     * @param x     value to be assigned
     */
    public final void setAdd( int i, int j, double x ) {
        a[r+i*ms+j*ns] += x;
    }

    /** Subtract a value to an entry of the matrix
     * @param i     row index (zero-based)
     * @param j     column index (zero-based)
     * @param x     value to be assigned
     */
    public final void setSub( int i, int j, double x ) {
        a[r+i*ms+j*ns] -= x;
    }

    /** Multiply a value to an entry of the matrix
     * @param i     row index (zero-based)
     * @param j     column index (zero-based)
     * @param x     value to be assigned
     */
    public final void setMul( int i, int j, double x ) {
        a[r+i*ms+j*ns] *= x;
    }

    /** Divide a value to an entry of the matrix
     * @param i     row index (zero-based)
     * @param j     column index (zero-based)
     * @param x     value to be assigned
     */
    public final void setDiv( int i, int j, double x ) {
        a[r+i*ms+j*ns] /= x;
    }

    /** swap two rows 
     * @param i1   row index 1 ( zero-based )
     * @param i2   row index 2 ( zero-based )
     */
    public final void swapRows( int i1, int i2 ) {
        double x;
        for ( int j=0; j<n; j++ )
        {
            x = get( i1, j );
            set( i1, j, get( i2, j ) );
            set( i2, j, x );
        }
    }

    /** swap two columns
     * @param i1   column index 1 ( zero-based )
     * @param i2   column index 2 ( zero-based )
     */
    public final void swapColumns( int j1, int j2 ) {
        double x;
        for ( int i=0; i<m; i++ )
        {
            x = get( i, j1 );
            set( i, j1, get( i, j2 ) );
            set( i, j2, x );
        }
    }

    /** Generate a copy of this matrix without share of data
     * @return      copy of this matrix
     */
    public final Matrix copy() {
        Matrix x = new Matrix( m, n );
        int c = 0, p = r, q = p;

        for ( int j=0; j<n; j++ )
        {
            q = p;
            for ( int i=0; i<m; i++ )
            {
                x.a[c++] = a[q];
                q += ms;
            }
            p += ns;
        }

        return x;
    }

    /** Implement the method in interface Cloneable
     * @return      copy of this matrix as an object
     */
    public Object clone() {
        return copy();
    }

    /** Let this matrix be as same as another matrix B (shallow copy)
     * @param b     matrix B
     * @return      this matrix
     */
    public Matrix refer( Matrix b ) {
        m = b.m;
        n = b.n;
        r = b.r;
        ms = b.ms;
        ns = b.ns;
        a = b.a;
        return this;
    }

    /** Copy all elements of matrix B to this matrix (deep copy)
     * @param b     Matrix B
     * @return      this matrix
     */
    public Matrix assign( Matrix b ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                set( i, j, b.get( i, j ) );
        return this;
    }

    /** Copy elements of matrix B to this matrix with bitmasks
     * @param bmk   Bitmask indicating whether permission of assignment
     * @param b     Matrix B
     */
    public Matrix assign( BitArray bmk, Matrix b ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ ) 
                if ( getBitMask( bmk, i, j ) )
                    set( i, j, b.get( i, j ) );
        return this;
    }

    /** Set all elements to a scalar
     * @param f     scalar
     * @return      this matrix
     */
    public Matrix assign( double f ) {
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                set( i, j, f );
        return this;
    }

    /** Copy elements to a scalar to this matrix with bitmasks
     * @param bmk   Bitmask indicating whether permission of assignment
     * @param b     Matrix B
     */
    public Matrix assign( BitArray bmk, double f ) {
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ ) 
                if ( getBitMask( bmk, i, j ) )
                    set( i, j, f );
        return this;
    }

    /** Get the submatrix of this matrix, sharing the same data
     * @param mr    row range
     * @param nr    column range
     * @return      sub-matrix
     */
    public Matrix slice( Range mr, Range nr ) {
        return new Matrix( this, mr, nr );

        /* This piece of code is no longer necessary

        Matrix x = new Matrix( this, mr, nr );

        int c1 = x.r;
        int c2 = x.r + (x.m-1)*x.ms;
        int c3 = x.r + (x.n-1)*x.ns;
        int c4 = x.r + (x.m-1)*x.ms + (x.n-1)*x.ns;

        if ( c1 < 0 || c1 >= x.a.length 
             || c2 < 0 || c2 >= x.a.length 
             || c3 < 0 || c3 >= x.a.length 
             || c4 < 0 || c4 >= x.a.length )
        {
            throw new IllegalArgumentException( 
                "Invalid range: "
                + m + "*" + n + " [" + mr.toString() + "," + nr.toString()
                + "]" );
        }
        return x;
        */
    }

    /* Order statistics: minimum
     * @return      the minimum value of a matrix
     */
    public double min() {
        double d = Double.MAX_VALUE;

        for ( int j=0; j<n; j++ ) 
            for ( int i=0; i<m; i++ )
            {
                if ( get(i,j) < d )
                    d = get(i,j);
            }
        return d;
    }

    /* Order statistics: maximum
     * @return      the maximum value of a matrix
     */
    public double max() {
        double d = Double.MIN_VALUE;

        for ( int j=0; j<n; j++ ) 
            for ( int i=0; i<m; i++ )
            {
                if ( get(i,j) > d )
                    d = get(i,j);
            }
        return d;
    }

    /* Order statistics: minimum of a column
     * @param column the column number
     * @return      the minimum value of the matrix column
     */
    public double min( int column ) {
        double d = Double.MAX_VALUE;

        for ( int i=0; i<m; i++ )
        {
            if ( get(i,column) < d )
                d = get(i,column);
        }
        return d;
    }

    /* Order statistics: maximum of a column
     * @param column the column number
     * @return      the maximum value of the matrix column
     */
    public double max( int column ) {
        double d = Double.MIN_VALUE;

        for ( int i=0; i<m; i++ )
        {
            if ( get(i,column) > d )
                d = get(i,column);
        }
        return d;
    }

    /** Transpose this matrix, sharing the same data
     * @return      new tranposed matrix
     */
    public Matrix transpose() {
        Matrix x = new Matrix( this );
        int t;
        t = x.m; x.m = x.n; x.n = t;
        t = x.ms; x.ms = x.ns; x.ns = t;
        return x;
    }

    /** Flip vertically this matrix, sharing the same data
     * @return      new flipped matrix
     */
    public Matrix flip() {
        Matrix x = new Matrix( this );
        x.r += (x.m - 1) * x.ms;
        x.ms = -x.ms;
        return x;
    }

    /** Mirror horizontally this matrix, sharing the same data
     * @return      new mirrored matrix
     */
    public Matrix mirror() {
        Matrix x = new Matrix( this );
        x.r += (x.n - 1) * x.ns;
        x.ns = -x.ns;
        return x;
    }

    /** Negate all elements of a matrix
     * @return      this matrix
     */
    public Matrix selfneg() {
        for ( int j=0; j<n; j++ ) 
            for ( int i=0; i<m; i++ )
                set( i, j, -get(i,j) );
        return this;
    }

    /** Generate a new matrix that is the negative of this matrix
     * @return      new negative matrix
     */
    public Matrix uminus() {
        Matrix x = new Matrix( m, n );
        for ( int j=0; j<n; j++ ) 
            for ( int i=0; i<m; i++ )
                x.set( i, j, -get(i,j) );
        return x;
    }
    
    /** Compute the sum of this matrix and another matrix (C=A+B)
     * @param b     second operand B
     * @return      new sum matrix C
     */
    public Matrix plus( Matrix b ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        Matrix x = new Matrix( m, n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                x.set( i, j, get(i,j) + b.get(i,j) );
        return x;
    }

    /** Add another matrix to this matrix (A+=B)
     * @param b     matrix B
     * @return      this matrix, which is the sum
     */
    public Matrix selfadd( Matrix b ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                setAdd( i, j, b.get(i,j) );
        return this;
    }

    /** Compute the difference between this matrix and another matrix (C=A-B)
     * @param b     matrix B
     * @return      new difference matrix C
     */
    public Matrix minus( Matrix b ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        Matrix x = new Matrix( m, n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                x.set( i, j, get(i,j) - b.get(i,j) );
        return x;
    }

    /** Subtract another matrix from this matrix (A-=B)
     * @param b     matrix B
     * @return      this matrix, which is the difference
     */
    public Matrix selfsub( Matrix b ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                setSub( i, j, b.get(i,j) );
        return this;
    }

    /** Compute the product of this matrix by a scalar (C=fA)
     * @param f     scalar operand
     * @return      new product matrix 
     */
    public Matrix times( double f ) {
        Matrix x = new Matrix( m, n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                x.set( i, j, get(i,j) * f );
        return x;
    }

    /** Compute elemental product of this matrix and another matrix (C=A.*B)
     * @param b     second operand matrix B
     * @return      new product matrix
     */
    public Matrix etimes( Matrix b ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        Matrix x = new Matrix( m, n ); 
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                x.set( i, j, get(i,j) * b.get(i,j) );
        return x;
    }

    /** Compute the product of this matrix and another matrix (C=AB)
     * @param b     second operand matrix B
     * @return      new product matrix
     */
    public Matrix times( Matrix b ) {
        if ( n != b.m )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        Matrix x = new Matrix( m, b.n );
        for ( int j=0; j<b.n; j++ )
            for ( int i=0; i<m; i++ )
            {
                double t = 0.0;
                for ( int k=0; k<n; k++ )
                    t += get( i, k ) * b.get( k, j );
                x.set( i, j, t );
            }
        return x;
    }

    /** Multiply this matrix by a scalar (A*=f)
     * @param f     scalar 
     * @return      this matrix
     */
    public Matrix selfmul( double f ) {
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                setMul( i, j, f );
        return this;
    }

    /** Elementally multiply this matrix by another matrix (A.*=B)
     * @param b     matrix B
     * @return      this matrix
     */
    public Matrix selfemul( Matrix b ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                setMul( i, j, b.get(i,j) );
        return this;
    }

    /** Left multiply by another matrix (at right) (A=BA)
     * @param b     matrix B
     * @return      this matrix
     */
    public Matrix selflmul( Matrix b ) {
        return refer( b.times( this ) );
    }

    /** Right multiply  by another matrix (at left) (A=AB)
     * @param b     matrix B
     * @return      this matrix
     */
    public Matrix selfrmul( Matrix b ) {
        return refer( times( b ) );
    }

    /** Compute elemental division of this matrix and another matrix (C=A./B)
     * @param b     matrix B
     * @return      new quotient matrix
     */
    public Matrix edivide( Matrix b ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        Matrix x = new Matrix( m, n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                x.set( i, j, get(i,j) / b.get(i,j) );
        return x;
    }

    /** Compute left division of this matrix and another matrix (C=B^-1A)
     * @param b     matrix B
     * @return      new quotient matrix
     */
    public Matrix ldivide( Matrix b ) {
        return b.solve( this );
    }

    /** Compute left division of this matrix and another matrix (C=AB^-1)
     * @param b     matrix B
     * @return      new quotient matrix
     */
    public Matrix rdivide( Matrix b ) {
        return b.transpose().solve( transpose() ).transpose();
    }

    /** Elementally divide this matrix by another matrix (A=A./B)
     * @param b     matrix B
     * @return      this matrix
     */
    public Matrix selfediv( Matrix b ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                setDiv( i, j, b.get(i,j) );
        return this;
    }

    /** Left divide another matrix B from this matrix (A=B^-1A)
     * @param b     matrix B
     * @return      new quotient matrix
     */
    public Matrix selfldiv( Matrix b ) {
        return refer( ldivide( b ) );
    }

    /** Right divide another matrix B from this matrix (A=AB^-1)
     * @param b     matrix B
     * @return      new quotient matrix
     */
    public Matrix selfrdiv( Matrix b ) {
        return refer( rdivide( b ) );
    }

    /** Return a bit array indicating whether elements of this matrix is
     * greater than a scalar value
     * @param f     a scalar
     * @return      new bit array
     */
    public BitArray gt( double f ) {
        BitArray bmk = new BitArray( m*n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                if ( get(i,j) > f )
                    setBitMask( bmk, i, j );
        return bmk;
    }

    /** Return a bit array indicating whether elements of this matrix is
     * greater than corresponding elements of another matrix B
     * @param b     matrix B
     * @return      new bit array
     */
    public BitArray gt( Matrix b ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        BitArray bmk = new BitArray( m*n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                if ( get(i,j) > b.get(i,j) )
                    setBitMask( bmk, i, j );
        return bmk;
    }

    /** Return a bit array indicating whether elements of this matrix is
     * greater than or equal to a scalar 
     * @param f     a scalar
     * @return      new bit array
     */
    public BitArray ge( double f ) {
        BitArray bmk = new BitArray( m*n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                if ( get(i,j) >= f )
                    setBitMask( bmk, i, j );
        return bmk;
    }

    /** Return a bit array indicating whether elements of this matrix is
     * greater than or equal to corresponding elements of another matrix B
     * @param b     matrix B
     * @return      new bit array
     */
    public BitArray ge( Matrix b ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        BitArray bmk = new BitArray( m*n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                if ( get(i,j) >= b.get(i,j) )
                    setBitMask( bmk, i, j );
        return bmk;
    }


    /** Return a bit array indicating whether elements of this matrix is
     * less than a scalar
     * @param f     a scalar
     * @return      new bit array
     */
    public BitArray lt( double f ) {
        BitArray bmk = new BitArray( m*n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                if ( get(i,j) < f )
                    setBitMask( bmk, i, j );
        return bmk;
    }

    /** Return a bit array indicating whether elements of this matrix is
     * less than corresponding elements of another matrix B
     * @param b     matrix B
     * @return      new bit array
     */
    public BitArray lt( Matrix b ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        BitArray bmk = new BitArray( m*n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                if ( get(i,j) < b.get(i,j) )
                    setBitMask( bmk, i, j );
        return bmk;
    }

    /** Return a bit array indicating whether elements of this matrix is
     * less than or equal to a scalar
     * @param f     a scalar
     * @return      new bit array
     */
    public BitArray le( double f ) {
        BitArray bmk = new BitArray( m*n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                if ( get(i,j) <= f )
                    setBitMask( bmk, i, j );
        return bmk;
    }

    /** Return a bit array indicating whether elements of this matrix is
     * less than or equal to corresponding elements of another matrix B
     * @param b     matrix B
     * @return      new bit array
     */
    public BitArray le( Matrix b ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        BitArray bmk = new BitArray( m*n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                if ( get(i,j) <= b.get(i,j) )
                    setBitMask( bmk, i, j );
        return bmk;
    }

    /** Return a bit array indicating whether elements of this matrix is
     * equal to a scalar
     * @param f     a scalar
     * @return      new bit array
     */
    public BitArray eq( double f ) {
        BitArray bmk = new BitArray( m*n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                if ( get(i,j) == f )
                    setBitMask( bmk, i, j );
        return bmk;
    }

    /** Return a bit array indicating whether elements of this matrix is
     * equal to corresponding elements of another matrix B
     * @param b     matrix B
     * @return      new bit array
     */
    public BitArray eq( Matrix b ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        BitArray bmk = new BitArray( m*n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                if ( get(i,j) == b.get(i,j) )
                    setBitMask( bmk, i, j );
        return bmk;
    }

    /** Return a bit array indicating whether elements of this matrix is
     * equal to a scalar within a given error
     * @param f     a scalar
     * @param eps   error below which elements can be considered equal.
     * @return      new bit array
     */
    public BitArray eq( double f, double eps ) {
        BitArray bmk = new BitArray( m*n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                if ( Math.abs( get(i,j) - f ) < eps )
                    setBitMask( bmk, i, j );
        return bmk;
    }

    /** Return a bit array indicating whether elements of this matrix is
     * equal to corresponding elements of another matrix B within a given error
     * @param b     matrix B
     * @param eps   error below which elements can be considered equal.
     * @return      new bit array
     */
    public BitArray eq( Matrix b, double eps ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        BitArray bmk = new BitArray( m*n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                if ( Math.abs( get(i,j) - b.get(i,j) ) < eps )
                    setBitMask( bmk, i, j );
        return bmk;
    }

    /** Return a bit array indicating whether elements of this matrix is not
     * equal to a scalar
     * @param f     a scalar
     * @return      new bit array
     */
    public BitArray ne( double f ) {
        BitArray bmk = new BitArray( m*n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                if ( get(i,j) != f )
                    setBitMask( bmk, i, j );
        return bmk;
    }

    /** Return a bit array indicating whether elements of this matrix is not
     * equal to corresponding elements of another matrix B
     * @param b     matrix B
     * @return      new bit array
     */
    public BitArray ne( Matrix b ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        BitArray bmk = new BitArray( m*n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                if ( get(i,j) != b.get(i,j) )
                    setBitMask( bmk, i, j );
        return bmk;
    }

    /** Return a bit array indicating whether elements of this matrix is not
     * equal to a scalar beyond a given error
     * @param f     a scalar
     * @param eps   error below which elements can be considered equal.
     * @return      new bit array
     */
    public BitArray ne( double f, double eps ) {
        BitArray bmk = new BitArray( m*n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                if ( Math.abs( get(i,j) - f ) >= eps )
                    setBitMask( bmk, i, j );
        return bmk;
    }

    /** Return a bit array indicating whether elements of this matrix is not
     * equal to corresponding elements of another matrix B beyond a given error
     * @param b     matrix B
     * @param eps   error below which elements can be considered equal.
     * @return      new bit array
     */
    public BitArray ne( Matrix b, double eps ) {
        if ( m != b.m || n != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        BitArray bmk = new BitArray( m*n );
        for ( int j=0; j<n; j++ )
            for ( int i=0; i<m; i++ )
                if ( Math.abs( get(i,j) - b.get(i,j) ) >= eps )
                    setBitMask( bmk, i, j );
        return bmk;
    }

    /** Solve the linear equations that AX=B. 
     * @param b     matrix B
     * @return      solution X
     */
    public Matrix solve( Matrix b ) {
        if ( b.m != m )
            throw new IllegalArgumentException( dimErrMsg( b ) );
        Matrix x = new Matrix( b.m, b.n );
        if ( m == n )
        {
            int [] index = new int[m];
            Matrix lum = copy();
            lum.lud( index );

            // exchange rows
            for ( int k=0; k<n; k++ )
            {
                for ( int j=0; j<b.n; j++ )
                    x.set( k, j, b.get(index[k],j) );
            }

            // solve L*X=B
            for ( int k=0; k<n; k++ )
            {
                for ( int i=k+1; i<n; i++ )
                    for ( int j=0; j<b.n; j++ )
                        x.setSub( i, j, 
                                  x.get(k,j)*lum.get(i,k) );
            }

            // solve U*X=X
            for ( int k=n-1; k>=0; k-- ) 
            {
                double d = lum.get( k, k );
                if ( d == 0 )
                    throw new RuntimeException( "Singular" );
                d = 1.0 / d;
                for ( int j=0; j<b.n; j++ )
                    x.setMul( k, j, d );
                for ( int i=0; i<k; i++ )
                    for ( int j=0; j<b.n; j++ )
                        x.setSub( i, j, 
                                  x.get(k,j)*lum.get(i,k) );
            }
        }
        else
            return null;
        return x;
    }

    /** compute the determinant of a matrix
     * @return      the determinant of this matrix
     */
    public double det() {
        if ( m != n )
            throw new IllegalArgumentException( dimErrMsg() );
        int [] index = new int[m];
        Matrix lum = copy();
        double x = (double) lum.lud( index );
        for ( int j=0; j<n; j++ )
            x *= lum.get( j, j );
        return x;
    }

    public int rank() {
        throw new UnsupportedOperationException( "rank" );
    }

    public int rank( double eps ) {
        throw new UnsupportedOperationException( "rank" );
    }

    /** invert a matrix
     * @return    the inverse of this matrix
     */
    public Matrix invert() {
        return solve( eye( n ) );
    }

    /** LU decompose this matrix using Crout's algorithm.  Note that
     * L and U are saved in this matrix as output.
     * @param index   output: original columns' corresponding indexes 
     * in the decomposed LU matrix
     * @return        +1 or -1 indicating the sign for det()
     */
    public int lud( int[] index ) {
        if ( m < n )
            throw new IllegalArgumentException( dimErrMsg() );

        int sign = 1;
        for ( int i=0; i<m; i++ )
            index[i] = i;

        // for each column do
        for ( int j=0; j<n; j++ ) 
        {
            double x, sum;
            for ( int i=0; i<j; i++ )
            {
                sum = get( i, j );
                for ( int k=0; k<i; k++ )
                    sum -= get( i, k ) * get( k, j );
                set( i, j, sum );
            }

            double maxv = 0.0;
            int maxidx = -1;

            for ( int i=j; i<m; i++ )
            {
                sum = get( i, j );
                for ( int k=0; k<j; k++ )
                    sum -= get( i, k ) * get( k, j );
                set( i, j, sum );
                x = Math.abs( sum );
                if ( x >= maxv )
                {
                    maxv = x;
                    maxidx = i;
                }
            }
            
            if ( maxidx != j )
            {
                int t = index[j];
                index[j] = index[maxidx];
                index[maxidx] = t;

                sign = -sign;
                swapRows( j, maxidx );
            }

            if ( get( j, j ) != 0.0 ) 
            {
                x = 1.0 / get( j, j );
                for ( int i=j+1; i<n; i++ )
                    setMul( i, j, x );
            }
        }

        return sign;
    }

    /** Strassen's matrix multiplication algorithm (test purpose only)
     * @param b     the second matrix B
     * @return      product of this matrix and matrix B
     */
    public final Matrix strassen( Matrix b ) {
        if ( m != n || n != b.m || b.m != b.n )
            throw new IllegalArgumentException( dimErrMsg( b ) );

        if ( 1 == n )
            return new Matrix( 1, 1, get(0,0) * b.get(0,0) );

        if ( 0 != (n&1) )
            throw new IllegalArgumentException( 
                "Strassen's algorithm is not applicable for matrix "
                + n + "*" + n );

        Range r1 = new Range( 0, n>>1, 1 );
        Range r2 = new Range( n>>1, n>>1, 1 );

        Matrix x = new Matrix( n, n );
            
        Matrix a11 = slice( r1, r1 );
        Matrix a12 = slice( r1, r2 );
        Matrix a21 = slice( r2, r1 );
        Matrix a22 = slice( r2, r2 );
    
        Matrix b11 = b.slice( r1, r1 );
        Matrix b12 = b.slice( r1, r2 );
        Matrix b21 = b.slice( r2, r1 );
        Matrix b22 = b.slice( r2, r2 );

        /* p1 = (a11+a22)*(b11+b22) */
        Matrix p1 = a11.plus( a22 ).strassen( b11.plus( b22 ) );
        
        /* p2 = (a21+a22)*b11 */
        Matrix p2 = a21.plus( a22 ).strassen( b11 );
        
        /* p3 = a11*(b12-b22) */
        Matrix p3 = a11.strassen( b12.minus( b22 ) );
        
        /* p4 = a22*(-b11+b21) */
        Matrix p4 = a22.strassen( b21.minus( b11 ) );
        
        /* p5 = (a11+a12)*b22 */
        Matrix p5 = a11.plus( a12 ).strassen( b22 );
        
        /* p6 = (-a11+a21)*(b11+b12) */
        Matrix p6 = a21.minus( a11 ).strassen( b11.plus( b12 ) );
        
        /* p7 = (a12-a22)*(b21+b22) */
        Matrix p7 = a12.minus( a22 ).strassen( b21.plus( b22 ) );
        
        /* c11 = p1+p4-p5+p7 */
        x.slice(r1,r1).assign(p1).selfadd(p4).selfsub(p5).selfadd(p7);
    
        /* c12 = p3+p5 */
        x.slice(r1,r2).assign(p3).selfadd(p5);
    
        /* c21 = p2+p4 */
        x.slice(r2,r1).assign(p2).selfadd(p4);
    
        /* c22 = p1+p3-p2+p6 */
        x.slice(r2,r2).assign(p1).selfadd(p3).selfsub(p2).selfadd(p6);

        return x;
    }

    /** binary file of unsigned bytes */
    public final static int DT_BYTE = 1;
    /** binary file of signed short integers (16-bits) */
    public final static int DT_SHORT = 2;
    /** binary file of signed integers (32-bits) */
    public final static int DT_INT = 3;
    /** binary file in IEEE floating format (32-bits) */
    public final static int DT_FLOAT = 4;
    /** binary file in IEEE floating format (64-bits) */
    public final static int DT_DOUBLE = 5;

    private static int[] dt_size = new int[]{ 0, 1, 2, 4, 4, 8 };

    /** read a matrix from a binary file.  Note at least one dimension should
     * be positive
     * @param fn      the file name
     * @param type    data type: DT_BYTE, DT_SHORT, DT_INT, DT_FLOAT, DT_DOUBLE
     * @param m       number of rows (determined by file size if <= 0)
     * @param n       number of columns (determined by file size if <=0)
     * @param skip    number of bytes to be skipped from the beginning
     */
    public static Matrix load( String fn, int type, int m, int n, long skip ) 
        throws IOException {

        if ( m <= 0 && n <= 0 )
            return null;

        RandomAccessFile file = new RandomAccessFile( fn, "r" );
        long l = file.length() - skip;
        if ( l < 0 )
        {
            file.close();
            return null;
        }

        l /= dt_size[type];
        if ( m <= 0 )
            m = (int)( l / n );
        else if ( n <= 0 )
            n = (int)( l / m );
        if ( m <= 0 || n <= 0 )
        {
            file.close();
            return null;
        }
        
        if ( skip != 0 )
            file.seek( skip );

        Matrix x = new Matrix( m, n );
        
        /* file is default saved by columns */
        switch ( type )
        {
        case DT_BYTE:  /* Java byte is signed!! */
            for ( int j=0; j<n; j++ )
                for ( int i=0; i<m; i++ )
                    x.set( i, j, (double) (0xff & (int) file.readByte()) );
            break;
        case DT_SHORT:
            for ( int j=0; j<n; j++ )
                for ( int i=0; i<m; i++ )
                    x.set( i, j, (double) file.readShort() );
            break;
        case DT_INT:
            for ( int j=0; j<n; j++ )
                for ( int i=0; i<m; i++ )
                    x.set( i, j, (double) file.readInt() );
            break;
        case DT_FLOAT:
            for ( int j=0; j<n; j++ )
                for ( int i=0; i<m; i++ )
                    x.set( i, j, (double) file.readFloat() );
            break;
        case DT_DOUBLE:
            for ( int j=0; j<n; j++ )
                for ( int i=0; i<m; i++ )
                    x.set( i, j, file.readDouble() );
            break;
        default:
            file.close();
            throw new IllegalArgumentException( "Illegal data type" );
        }
        file.close();

        return x;
    }

    /** read a matrix from a binary file.  Note at least one dimension should
     * be positive
     * @param fn      the file name
     * @param type    data type: DT_BYTE, DT_SHORT, DT_INT, DT_FLOAT, DT_DOUBLE
     * @param m       number of rows (determined by file size if <= 0)
     * @param n       number of columns (determined by file size if <=0)
     */
    public static Matrix load( String fn, int type, int m, int n )
        throws IOException {
        return load( fn, type, m, n, 0 );
    }

    /** write a matrix to a binary file.  
     * @param fn      the file name
     * @param type    data type: DT_BYTE, DT_SHORT, DT_INT, DT_FLOAT, DT_DOUBLE
     * @param skip    number of bytes to be skipped from the beginning
     */
    public void save( String fn, int type, long skip ) 
        throws IOException {

        RandomAccessFile file = new RandomAccessFile( fn, "rw" );
        if ( skip != 0 )
            file.seek( skip );
            
        /* file is default saved by columns */
        switch ( type )
        {
        case DT_BYTE:
            for ( int j=0; j<n; j++ )
                for ( int i=0; i<m; i++ )
                    file.writeByte( (byte) (int) get( i, j ) );
            break;
        case DT_SHORT:
            for ( int j=0; j<n; j++ )
                for ( int i=0; i<m; i++ )
                    file.writeShort( (short) (int) get( i, j ) );
            break;
        case DT_INT:
            for ( int j=0; j<n; j++ )
                for ( int i=0; i<m; i++ )
                    file.writeInt( (int) get( i, j ) );
            break;
        case DT_FLOAT:
            for ( int j=0; j<n; j++ )
                for ( int i=0; i<m; i++ )
                    file.writeFloat( (float) get( i, j ) );
            break;
        case DT_DOUBLE:
            for ( int j=0; j<n; j++ )
                for ( int i=0; i<m; i++ )
                    file.writeDouble( get( i, j ) );
            break;
        default:
            file.close();
            throw new IllegalArgumentException( "Illegal data type" );
        }
        file.close();
    }

    /** write a matrix to a binary file.  
     * @param fn      the file name
     * @param type    data type: DT_BYTE, DT_SHORT, DT_INT, DT_FLOAT, DT_DOUBLE
     */
    public void save( String fn, int type ) throws IOException {
        save( fn, type, 0 );
    }

    /** print a matrix with specified format to the standard output
     * @param width    field width of each elements
     * @param pecision number of fraction digits
     * @param indent   number of additional spaces at the beginning of a line
     */
    public void print( int width, int precision, int indent ) {
        print( new PrintWriter( System.out, true ), width, precision, indent );
    }

    /** print a matrix with user-defined number format to the standard output
     * @param format   user-defined number format
     * @param width    field width of each elements
     * @param indent   number of additional spaces at the beginning of a line
     */
    public void print( NumberFormat format, int width, int indent ) {
        print( new PrintWriter( System.out, true ), format, width, indent );
    } 

    /** print a matrix with specified format
     * @param output   the printing device
     * @param width    field width of each elements
     * @param pecision number of fraction digits
     * @param indent   number of additional spaces at the beginning of a line
     */
    public void print( PrintWriter output, int width, int precision, 
                       int indent ) {
        DecimalFormat fmt = new DecimalFormat();
        fmt.setDecimalFormatSymbols( new DecimalFormatSymbols(Locale.US) );
        fmt.setMinimumIntegerDigits( 1 );
        fmt.setMaximumFractionDigits( precision );
        fmt.setMinimumFractionDigits( precision );
        fmt.setGroupingUsed( false );
        print( output, fmt, width, indent );
    }

    /** print a matrix with user-defined number format
     * @param output   the printing device
     * @param format   user-defined number format
     * @param width    field width of each elements
     * @param indent   number of additional spaces at the beginning of a line
     */
    public void print( PrintWriter output, NumberFormat format, 
                       int width, int indent ) {
        for ( int i=0; i<m; i++ )
        {
            for ( int j=0; j<n; j++ )
            {
                String str = format.format( get(i,j) );
                int space = width - str.length();
                if ( 0 == j )
                {
                    space += indent;
                }
                else
                {
                    output.print( ' ' );
                    space--;
                }
                for ( ; space>0; space-- )
                    output.print( ' ' );
                output.print( str );
            }
            output.println();
        }
    }
}

