
public final class VecMath
{
    private VecMath() {}


    public static void vpv(double result[], double v0[], double v1[])
    {
        for (int i = (v0.length)-1; (i) >= 0; --i)
            result[i] = v0[i] + v1[i];
    }
    public static void vpv(int result[], int v0[], int v1[])
    {
        for (int i = (v0.length)-1; (i) >= 0; --i)
            result[i] = v0[i] + v1[i];
    }
    public static void vmv(double result[], double v0[], double v1[])
    {
        for (int i = (v0.length)-1; (i) >= 0; --i)
            result[i] = v0[i] - v1[i];
    }

    public static void vxs(double result[], double v[], double s)
    {
        for (int i = (v.length)-1; (i) >= 0; --i)
            result[i] = v[i] * s;
    }

    public static void mxs(double result[][], double M[][], double s)
    {
        for (int i = (M.length)-1; (i) >= 0; --i)
            if (result[i] != null)
                vxs(result[i], M[i], s);
    }
    public static void mxs(double result[][][], double M[][][], double s)
    {
        for (int i = (M.length)-1; (i) >= 0; --i)
            mxs(result[i], M[i], s);
    }
    public static void mxs(double result[][][][], double M[][][][], double s)
    {
        for (int i = (M.length)-1; (i) >= 0; --i)
            mxs(result[i], M[i], s);
    }

    public static void sxv(double result[], double s, double v[])
    {
        for (int i = (v.length)-1; (i) >= 0; --i)
            result[i] = s * v[i];
    }


    public static void xv2(double result[ ], double v[ ])
    {
        result[0] = -v[1];
        result[1] = v[0];
    }

    public static double vxv2(double v[ ], double w[ ])
    {
        return v[0]*w[1] - v[1]*w[0];
    }
    public static void vxv3(double v0[ ], final double v1[ ], final double v2[ ])
    {
        v0[0] = v1[1]*v2[2] - v1[2]*v2[1];
        v0[1] = v1[2]*v2[0] - v1[0]*v2[2];
        v0[2] = v1[0]*v2[1] - v1[1]*v2[0];
    }

    public static double vxvxv3(double v0[ ], double v1[ ], double v2[ ])
    {
        return v0[0] * (v1[1]*v2[2] - v1[2]*v2[1])
             + v0[1] * (v1[2]*v2[0] - v1[0]*v2[2])
             + v0[2] * (v1[0]*v2[1] - v1[1]*v2[0]);
    }
    public static void crossprod(double result[ ], double vectors[ ][ ])
    {
        int dim = result.length;
        { if (!(vectors.length == dim-1)) HyperbolicUtils.set_exception("VecMath.java","crossprod" , "vectors.length == dim-1" ); };
        if (vectors.length > 0)
            { if (!(vectors[0].length == dim)) HyperbolicUtils.set_exception("VecMath.java","crossprod" , "vectors[0].length == dim" ); };

        if (true)
        {

            switch(dim)
            {
                case 1:
                    result[0] = 1.;
                    return;
                case 2:
                    xv2(result, vectors[0]);
                    return;
                case 3:

                    vxv3(result, vectors[0], vectors[1]);
                    return;
            }
        }
        else
        {
            if (dim == 1 || dim == 2 || dim == 3)
                System.err.println("WARNING: doing slow crossprod for testing, dim = "+dim+"");
        }

        {
            double cofactorMatrix[][] = new double[dim-1][dim-1];
            int sign = ((dim-1) % 2 == 0 ? 1 : -1);
            for (int i = 0; (i) < (dim); ++i)
            {

                {
                    for (int j = 0; (j) < (dim-1); ++j)
                    for (int k = 0; (k) < (dim-1); ++k)
                        cofactorMatrix[j][k] = vectors[j][k>=i?k+1:k];
                }
                result[i] = detDestructive(cofactorMatrix) * sign;
                sign = -sign;
            }
        }
    }


    public static void vpsxv(double result[], double v0[],
                                              double s1, double v1[])
    {
        for (int i = (v0.length)-1; (i) >= 0; --i)
            result[i] = v0[i] + s1*v1[i];
    }

    public static void vpsxv(int result[], int v0[],
                                           int s1, int v1[])
    {
        for (int i = (v0.length)-1; (i) >= 0; --i)
            result[i] = v0[i] + s1*v1[i];
    }

    public static void vpvmv(double result[], double v0[], double v1[], double v2[])
    {
        for (int i = (v0.length)-1; (i) >= 0; --i)
            result[i] = v0[i] + (v1[i] - v2[i]);
    }

    public static void sxvpsxv(double result[], double s0, double v0[],
                                                double s1, double v1[])
    {
        for (int i = (v0.length)-1; (i) >= 0; --i)
            result[i] = s0*v0[i] + s1*v1[i];
    }


    public static void mpv(double result[][], double M[][], double v[])
    {
        for (int i = (M.length)-1; (i) >= 0; --i)
            if (M[i] != null)
                vpv(result[i], M[i], v);
    }
    public static void mpv(double result[][][], double M[][][], double v[])
    {
        for (int i = (M.length)-1; (i) >= 0; --i)
            mpv(result[i], M[i], v);
    }
    public static void mpv(double result[][][][], double M[][][][], double v[])
    {
        for (int i = (M.length)-1; (i) >= 0; --i)
            mpv(result[i], M[i], v);
    }

    public static void mmv(double result[][], double M[][], double v[])
    {
        for (int i = (M.length)-1; (i) >= 0; --i)
            vmv(result[i], M[i], v);
    }


    public static double sqrdTwiceTriangleArea(double a[],
                                               double b[],
                                               double c[])
    {
        double ab_dot_ab = 0.;
        double ab_dot_ac = 0.;
        double ac_dot_ac = 0.;
        for (int i = (a.length)-1; (i) >= 0; --i)
        {
            double ai = a[i];
            double ab = b[i] - ai;
            double ac = c[i] - ai;
            ab_dot_ab += ab*ab;
            ab_dot_ac += ab*ac;
            ac_dot_ac += ac*ac;
        }
        return ab_dot_ab*ac_dot_ac - ((ab_dot_ac)*(ab_dot_ac));
    }



    public static void lerp(double result[], double p0[], double p1[], double t)
    {
        for (int i = (p0.length)-1; (i) >= 0; --i)
            result[i] = p0[i] + t*(p1[i]-p0[i]);
    }
    public static void bary(double result[],
                            double p0[],
                            double p1[], double t1,
                            double p2[], double t2)
    {
        for (int i = (p0.length)-1; (i) >= 0; --i)
            result[i] = (1.-t1-t2)*p0[i] + t1*p1[i] + t2*p2[i];
    }

    public static void cerp(double result[],
                            double p0[], double p1[],
                            double v0[], double v1[],
                            double t)
    {

        double t2 = t*t;
        double coeff_p0 = (2*t - 3)*t2 + 1;
        double coeff_p1 = (-2*t + 3)*t2;
        double coeff_v0 = ((t - 2)*t + 1)*t;
        double coeff_v1 = (t - 1)*t2;

        if (false)
        {
            coeff_p0 = (1-t);
            coeff_p1 = t;
            coeff_v0 = 0;
            coeff_v1 = 0;
        }

        sxv(result, coeff_p0, p0);
        vpsxv(result, result, coeff_p1, p1);
        vpsxv(result, result, coeff_v0, v0);
        vpsxv(result, result, coeff_v1, v1);
    }





    public static void interp5(double result[],
                               double p0[], double p1[],
                               double v0[], double v1[],
                               double a0[], double a1[],
                               double t)
    {

        double t2 = t*t, t3 = t2*t, t4 = t3*t, t5 = t4*t;
        double tVec[] = {1,t,t2,t3,t4,t5};
        final double M[][] = {
                { 1, 0, 0, 0, 0 , 0},
                { 0, 0, 1, 0, 0 , 0},
                { 0, 0, 0, 0, .5 , 0},
                {-10, 10, -6, -4, -1.5, .5},
                { 15, -15, 8, 7, 1.5, -1},
                { -6, 6, -3, -3, -.5 , .5},
        };
        double pva[][] = {p0,p1,v0,v1,a0,a1};


        double coeffs[] = vxm(tVec, M);
        vxm(result, coeffs, pva);

    }



    public static void copyvec(boolean result[], boolean v[])
    {
        for (int i = (((result.length)<=(v.length)?(result.length):(v.length)))-1; (i) >= 0; --i)
            result[i] = v[i];
    }

    public static void copyvec(int result[], int v[])
    {
        for (int i = (((result.length)<=(v.length)?(result.length):(v.length)))-1; (i) >= 0; --i)
            result[i] = v[i];
    }

    public static void copyvec(double result[], double v[])
    {
        for (int i = (((result.length)<=(v.length)?(result.length):(v.length)))-1; (i) >= 0; --i)
            result[i] = v[i];
    }

    public static void setcolumn(double M[][], int iCol, double v[])
    {
        for (int iRow = (((M.length)<=(v.length)?(M.length):(v.length)))-1; (iRow) >= 0; --iRow)
            M[iRow][iCol] = v[iRow];
    }

    public static void getcolumn(double v[], double M[][], int iCol)
    {
        for (int iRow = (((v.length)<=(M.length)?(v.length):(M.length)))-1; (iRow) >= 0; --iRow)
            v[iRow] = M[iRow][iCol];
    }

    public static void fillvec(double result[], double s)
    {
        for (int i = (result.length)-1; (i) >= 0; --i)
            result[i] = s;
    }

    public static void fillvec(int result[], int s)
    {
        for (int i = (result.length)-1; (i) >= 0; --i)
            result[i] = s;
    }

    public static void fillvec(boolean result[], boolean s)
    {
        for (int i = (result.length)-1; (i) >= 0; --i)
            result[i] = s;
    }

    public static void fillmat(double result[][], double s)
    {
        for (int i = (result.length)-1; (i) >= 0; --i)
            fillvec(result[i], s);
    }

    public static void fillmat(boolean result[][], boolean s)
    {
        for (int i = (result.length)-1; (i) >= 0; --i)
            fillvec(result[i], s);
    }

    public static void copymat(double result[][], double m[][])
    {
        for (int i = (((result.length)<=(m.length)?(result.length):(m.length)))-1; (i) >= 0; --i)
            for (int j = (((result[i].length)<=(m[i].length)?(result[i].length):(m[i].length)))-1; (j) >= 0; --j)
                result[i][j] = m[i][j];
    }

    public static void zerovec(double result[])
    {
        for (int i = (result.length)-1; (i) >= 0; --i)
            result[i] = (double)0;
    }

    public static void zerovec(int result[])
    {
        for (int i = (result.length)-1; (i) >= 0; --i)
            result[i] = (int)0;
    }
    public static void zeromat(double result[][])
    {
        for (int i = (result.length)-1; (i) >= 0; --i)
            zerovec(result[i]);
    }

    public static double dot(double v0[], double v1[])
    {
        double result = 0.;
        for (int i = (v0.length)-1; (i) >= 0; --i)
            result += v0[i]*v1[i];
        return result;
    }


    public static void vxm(double result[], double v[], double m[][])
    {
        { if (!(result != v)) HyperbolicUtils.set_exception("VecMath.java","vxm double" , "result != v" ); };
        zerovec(result);
        int i, vLen = v.length;
        for (i = 0; (i) < (vLen); ++i)
            vpsxv(result, result, v[i], m[i]);

        for (; i < m.length; ++i)
            vpv(result, result, m[i]);
    }


    public static void vxm(int result[], int v[], int m[][])
    {
        { if (!(result != v)) HyperbolicUtils.set_exception("VecMath.java","vxm int" , "result != v" ); };
        int vLen = v.length;
        zerovec(result);
        int i, j;
        for (i = 0; (i) < (vLen); ++i)
            vpsxv(result, result, v[i], m[i]);

        for (; i < m.length; ++i)
            vpv(result, result, m[i]);
    }

    public static void mxm(double result[ ][ ],
                           double m0[ ][ ],
                           double m1[ ][ ])
    {


        if (((m0).length==0 ? 0 : (m0)[0].length) == m1.length)
        {
            { if (!(result != m0 && result != m1)) HyperbolicUtils.set_exception("VecMath.java","mxm double" , "result != m0 && result != m1" ); };
            for (int i = (result.length)-1; (i) >= 0; --i)
                vxm(result[i], m0[i], m1);
        }
        else
        {


            if ((m0.length == m0[0].length || m0.length == m0[0].length+1)
             && m1.length == m0[0].length+1
             && m1[0].length == m0[0].length
             && result.length == m0[0].length+1
             && result[0].length == m0[0].length)
            {

                double M0[][] = identitymat(m1.length);
                for (int i = (m0.length)-1; (i) >= 0; --i)
                for (int j = (m0[i].length)-1; (j) >= 0; --j)
                    M0[i][j] = m0[i][j];
                mxm(result, M0, m1);
            }
            else
            {
             {System.out.print("    ");System.out.println("m0.length" + " = " + (m0.length));};
             {System.out.print("    ");System.out.println("((m0).length==0 ? 0 : (m0)[0].length)" + " = " + (((m0).length==0 ? 0 : (m0)[0].length)));};
             {System.out.print("    ");System.out.println("m1.length" + " = " + (m1.length));};
             {System.out.print("    ");System.out.println("((m1).length==0 ? 0 : (m1)[0].length)" + " = " + (((m1).length==0 ? 0 : (m1)[0].length)));};
             {System.out.print("    ");System.out.println("(m0.length == m0[0].length || m0.length == m0[0].length+1)" + " = " + ((m0.length == m0[0].length || m0.length == m0[0].length+1)));};
             {System.out.print("    ");System.out.println("m1.length == m0[0].length+1" + " = " + (m1.length == m0[0].length+1));};
             {System.out.print("    ");System.out.println("m1[0].length == m0[0].length" + " = " + (m1[0].length == m0[0].length));};
             {System.out.print("    ");System.out.println("result.length == m0[0].length+1" + " = " + (result.length == m0[0].length+1));};
             {System.out.print("    ");System.out.println("result[0].length == m0[0].length" + " = " + (result[0].length == m0[0].length));};
                { if (!(false)) HyperbolicUtils.set_exception("VecMath.java","mxm double" , "false" ); };
            }
        }

    }


    public static void mxm(int result[ ][ ],
                           int m0[ ][ ],
                           int m1[ ][ ])
    {
        { if (!(result != m0 && result != m1)) HyperbolicUtils.set_exception("VecMath.java","mxm int" , "result != m0 && result != m1" ); };
        for (int i = (result.length)-1; (i) >= 0; --i)
            vxm(result[i], m0[i], m1);
    }

    public static void mxv(double result[], double m[][], double v[])
    {
        { if (!(result != v)) HyperbolicUtils.set_exception("VecMath.java","mxv double" , "result != v" ); };
        int nRows = m.length;
        if (nRows == 0)
            return;
        int nCols = v.length;
        int i, j;

        for (i = 0; (i) < (nRows); ++i)
        {
            double m_i[] = m[i];
            double sum = 0.;
            for (j = 0; (j) < (nCols); ++j)
            {
                sum += m_i[j] * v[j];
            }
            result[i] = sum;
        }
    }


    public static void random(double v[])
    {
        do
            for (int i = (v.length)-1; (i) >= 0; --i)
                v[i] = (2.*Math.random() - 1.);
        while (normsqrd(v) > 1);
    }
    public static void random(double v[], java.util.Random generator)
    {
        do
            for (int i = (v.length)-1; (i) >= 0; --i)
                v[i] = (2.*generator.nextDouble() - 1.);
        while (normsqrd(v) > 1);
    }

    public static double distsqrd(double v0[], double v1[])
    {
        double result = 0.;
        for (int i = (v0.length)-1; (i) >= 0; --i)
        {
            double diff = v1[i]-v0[i];
            result += diff*diff;
        }
        return result;
    }
    public static double dist(double v0[], double v1[])
    {
        return Math.sqrt(distsqrd(v0, v1));
    }

    public static double normsqrd(double v[])
    {
        return dot(v,v);
    }
    public static double norm(double v[])
    {
        return Math.sqrt(dot(v,v));
    }

    public static int mini(double v[])
    {
        int mini = -1;
        double min = Double.POSITIVE_INFINITY;
        for (int i = (v.length)-1; (i) >= 0; --i)
        {
            if (v[i] < min)
            {
                min = v[i];
                mini = i;
            }
        }
        return mini;
    }
    public static int maxi(double v[])
    {
        int maxi = -1;
        double max = Double.NEGATIVE_INFINITY;
        for (int i = (v.length)-1; (i) >= 0; --i)
        {
            if (v[i] > max)
            {
                max = v[i];
                maxi = i;
            }
        }
        return maxi;
    }

    public static double min(double v[])
    {
        double min = Double.POSITIVE_INFINITY;
        for (int i = (v.length)-1; (i) >= 0; --i)
            min = Math.min(min, v[i]);
        return min;
    }
    public static double max(double v[])
    {
        double max = Double.NEGATIVE_INFINITY;
        for (int i = (v.length)-1; (i) >= 0; --i)
            max = Math.max(max, v[i]);
        return max;
    }

    public static double sum(double v[])
    {
        double sum = 0.;
        for (int i = (v.length)-1; (i) >= 0; --i)
            sum += v[i];
        return sum;
    }
    public static double average(double array[])
    {
        return sum(array) / array.length;
    }
    public static double geomAverage(double array[])
    {
        double logProd = 0.;
        for (int i = (array.length)-1; (i) >= 0; --i)
            logProd += Math.log(array[i]);
        return Math.exp(logProd / array.length);
    }

    public static void sum(double result[], double array[][])
    {
        zerovec(result);
        for (int i = (array.length)-1; (i) >= 0; --i)
            vpv(result, result, array[i]);
    }
    public static void average(double result[], double array[][])
    {
        sum(result, array);
        vxs(result, result, 1./array.length);
    }
    public static void sumIndexed(double result[], int inds[], double array[][])
    {
        zerovec(result);
        for (int i = (inds.length)-1; (i) >= 0; --i)
            vpv(result, result, array[inds[i]]);
    }
    public static void averageIndexed(double result[], int inds[], double array[][])
    {
        sumIndexed(result, inds, array);
        vxs(result, result, 1./inds.length);
    }
    public static void averageIndexed(double result[], int inds[][], double array[][])
    {
        int totalInds = 0;
        zerovec(result);
        for (int i = (inds.length)-1; (i) >= 0; --i)
        {
            int row[] = inds[i];
            for (int j = (row.length)-1; (j) >= 0; --j)
            {
                vpv(result, result, array[row[j]]);
                totalInds++;
            }
        }
        vxs(result, result, 1./totalInds);
    }


        private static void _bboxIndexed(double result[ ][],
                                         double coords[][],
                                         Object inds)
        {
            int nDims = coords[0].length;
            Class indsClass = inds.getClass();
            Class componentType = indsClass.getComponentType();
            { if (!(componentType != null)) HyperbolicUtils.set_exception("VecMath.java","_bboxIndexed" , "componentType != null" ); };
            if (componentType.isArray())
            {
                int n = java.lang.reflect.Array.getLength(inds);
                for (int i = (n)-1; (i) >= 0; --i)
                    _bboxIndexed(result,
                                 coords,
                                 java.lang.reflect.Array.get(inds, i));
            }
            else
            {
                int ints[] = (int[])inds;
                for (int iInt = (ints.length)-1; (iInt) >= 0; --iInt)
                for (int iDim = (nDims)-1; (iDim) >= 0; --iDim)
                {
                    double val = coords[ints[iInt]][iDim];
                    result[0][iDim] = Math.min(result[0][iDim], val);
                    result[1][iDim] = Math.max(result[1][iDim], val);
                }
            }
        }
    public static void bboxIndexed(double result[ ][],
                                   double coords[][],
                                   Object inds)
    {
        fillvec(result[0], Double.POSITIVE_INFINITY);
        fillvec(result[1], Double.NEGATIVE_INFINITY);
        _bboxIndexed(result, coords, inds);

    }

    public static void bbox(double result[ ][], double array[][])
    {
        for (int i = (result[0].length)-1; (i) >= 0; --i)
        {
            result[0][i] = Double.POSITIVE_INFINITY;
            result[1][i] = Double.NEGATIVE_INFINITY;
        }
        for (int i = (array.length)-1; (i) >= 0; --i)
        for (int j = (array[i].length)-1; (j) >= 0; --j)
        {
            result[0][j] = Math.min(result[0][j], array[i][j]);
            result[1][j] = Math.max(result[1][j], array[i][j]);
        }
    }

    public static void bboxUniform(double bbox[ ][], double array[][])
    {
        bbox(bbox, array);
        double bboxCenter[] = average(bbox);
        double bboxDims[] = vmv(bbox[1], bbox[0]);
        double bboxMaxDim = max(bboxDims);
        fillvec(bboxDims, bboxMaxDim);
        vpsxv(bbox[0], bboxCenter, -.5, bboxDims);
        vpsxv(bbox[1], bboxCenter, .5, bboxDims);
    }

    public static void bboxIntersect(double result[ ][],
                                     double bbox0[ ][],
                                     double bbox1[ ][])
    {
        int n = result[0].length;
        for (int i = 0; (i) < (n); ++i)
        {
            result[0][i] = ((bbox0[0][i])>=(bbox1[0][i])?(bbox0[0][i]):(bbox1[0][i]));
            result[1][i] = ((bbox0[1][i])<=(bbox1[1][i])?(bbox0[1][i]):(bbox1[1][i]));
        }
    }

    public static void bboxUnion(double result[ ][],
                                 double bbox0[ ][],
                                 double bbox1[ ][])
    {
        int n = result[0].length;
        for (int i = 0; (i) < (n); ++i)
        {
            result[0][i] = ((bbox0[0][i])<=(bbox1[0][i])?(bbox0[0][i]):(bbox1[0][i]));
            result[1][i] = ((bbox0[1][i])>=(bbox1[1][i])?(bbox0[1][i]):(bbox1[1][i]));
        }
    }

    public static boolean closedBBoxContains(double bbox[ ][],
                                             double point[],
                                             double eps)
    {
        for (int i = (point.length)-1; (i) >= 0; --i)
            if (!((bbox[0][i] - eps <=(point[i]))&&((point[i])<= bbox[1][i] + eps)))
                return false;
        return true;
    }
    public static boolean bboxInteriorContains(double bbox[ ][],
                                               double point[],
                                               double eps)
    {
        for (int i = (point.length)-1; (i) >= 0; --i)
            if (!((bbox[0][i] + eps <(point[i]))&&((point[i])< bbox[1][i] - eps)))
                return false;
        return true;
    }


    public static void transpose(double result[][], double M[][])
    {
        { if (!(result != M)) HyperbolicUtils.set_exception("VecMath.java","transpose" , "result != M" ); };
        for (int i = (M.length)-1; (i) >= 0; --i)
        for (int j = (M[i].length)-1; (j) >= 0; --j)
            result[j][i] = M[i][j];
    }

    public static int luDecompose(double A[ ][ ],
                             int n, int indx[])
        throws Exception
    {
        double vv[] = new double[n];
        int d = 1;
        for (int i = 0; (i) < (n); ++i)
        {

            double max = 0.;
            for (int j = 0; (j) < (n); ++j)
            {
                double temp = ((A[i][j]) < 0 ? -(A[i][j]) : (A[i][j]));
                if (temp > max)
                    max = temp;
            }
            if (max == 0.)
            {
                throw new Exception("Singular matrix in routine luDecompose");
            }
            vv[i] = 1./max;
        }
        for (int j = 0; (j) < (n); ++j)
        {
            for (int i = 0; (i) < (j); ++i)
            {
                double sum = A[i][j];
                for (int k = 0; (k) < (i); ++k)
                    sum -= A[i][k] * A[k][j];
                A[i][j] = sum;
            }

            double max = 0.;
            int imax = -1;
            for (int i = j; i < n; i++)
            {
                double sum = A[i][j];
                for (int k = 0; (k) < (j); ++k)
                    sum -= A[i][k] * A[k][j];
                A[i][j] = sum;
                double dum = vv[i] * ((sum) < 0 ? -(sum) : (sum));
                if (dum >= max)
                {
                    max = dum;
                    imax = i;
                }
            }
            { if (!(imax != -1)) HyperbolicUtils.set_exception("VecMath.java","luDecompose" , "imax != -1" ); };
            if (j != imax)
            {

                double temp[];
                {temp=(A[imax]);A[imax]=(A[j]);A[j]=(temp);};
                d = -d;
                vv[imax] = vv[j];
            }
            indx[j] = imax;
            if (A[j][j] == 0.)
            {





                A[j][j] = 1e-40;
            }
            if (j != n-1)
            {
                double invAjj = 1./A[j][j];
                for (int i = j+1; i < n; i++)
                    A[i][j] *= invAjj;
            }
        }
        return d;
    }

    public static void luBackSubstitute(double A[][],
                              int n,
                              int indx[],
                              double b[])
    {
        int ii = -1;






        for (int i = 0; (i) < (n); ++i)
        {
            int ip = indx[i];
            double sum = b[ip];
            b[ip] = b[i];
            if (ii >= 0)
                for (int j = ii; j < i; j++)
                    sum -= A[i][j] * b[j];
            else if (sum != 0.)
                ii = i;
            b[i] = sum;

        }
        for (int i = (n)-1; (i) >= 0; --i)
        {
            double sum = b[i];
            for (int j = i+1; j < n;j++)
                sum -= A[i][j] * b[j];
            b[i] = sum / A[i][i];
        }
    }

    public static void gaussj(double a[ ][ ], int n,
                              double b[ ][ ], int m)
        throws Exception
    {
        int indxc[], indxr[], ipiv[];
        int i,icol=-1,irow=-1,j,k,l,ll;
        double big,dum,pivinv,doubletemp;

        indxc = new int[n];
        indxr = new int[n];
        ipiv = new int[n];
        for (j = 0; (j) < (n); ++j)
            ipiv[j] = 0;
        for (i = 0; (i) < (n); ++i) {
            big = 0.;
            for (j = 0; (j) < (n); ++j) {
                if (ipiv[j] != 1) {
                    for (k = 0; (k) < (n); ++k) {
                        if (ipiv[k] == 0) {
                            if (((a[j][k]) < 0 ? -(a[j][k]) : (a[j][k])) >= big) {
                                big = ((a[j][k]) < 0 ? -(a[j][k]) : (a[j][k]));
                                irow = j;
                                icol = k;
                            }
                        } else if (ipiv[k] > 1) {

                            throw new Exception("Singular matrix in routine gaussj, 1");
                        }
                    }
                }
            }
            ++ipiv[icol];




            if (irow != icol) {
                for (l = 0; (l) < (n); ++l) {doubletemp=(a[irow][l]);a[irow][l]=(a[icol][l]);a[icol][l]=(doubletemp);};
                for (l = 0; (l) < (m); ++l) {doubletemp=(b[irow][l]);b[irow][l]=(b[icol][l]);b[icol][l]=(doubletemp);};
            }
            indxr[i]=irow;
            indxc[i]=icol;
            if (a[icol][icol] == 0.) {

                throw new Exception("Singular matrix in routine gaussj, 2");
            }
            pivinv = 1.0/a[icol][icol];
            a[icol][icol] = 1.;
            for (l = 0; (l) < (n); ++l) a[icol][l] *= pivinv;
            for (l = 0; (l) < (m); ++l) b[icol][l] *= pivinv;
            for (ll = 0; (ll) < (n); ++ll) {
                if (ll != icol) {
                    dum = a[ll][icol];
                    a[ll][icol] = 0.;
                    for (l = 0; (l) < (n); ++l) a[ll][l] -= a[icol][l]*dum;
                    for (l = 0; (l) < (m); ++l) b[ll][l] -= b[icol][l]*dum;
                }
            }
        }



        for (l = n-1; l >= 0; --l) {
            if (indxr[l] != indxc[l])
                for (k = 0; (k) < (n); ++k)
                    {doubletemp=(a[k][indxr[l]]);a[k][indxr[l]]=(a[k][indxc[l]]);a[k][indxc[l]]=(doubletemp);};
        }
    }

    public static void invertmat(double result[][], double M[][])
    {
        int nRows = M.length;
        int nCols = ((M).length==0 ? 0 : (M)[0].length);
        { if (!(nRows == result.length)) HyperbolicUtils.set_exception("VecMath.java","invertmat" , "nRows == result.length" ); };
        { if (!(nCols == ((result).length==0 ? 0 : (result)[0].length))) HyperbolicUtils.set_exception("VecMath.java","invertmat" , "nCols == ROWLENGTH(result)" ); };
        if (nRows == nCols)
        {




            double temp[][] = copymat(M);
            int indx[] = new int[nRows];
            try
            {
                          luDecompose(temp, nRows, indx);
            }
            catch (Exception e)
            {
                fillmat(result, Double.POSITIVE_INFINITY);
                return;
            }

            double col[] = new double[nRows];
            for (int j = 0; (j) < (nRows); ++j)
            {
                zerovec(col);
                col[j] = 1.;
                luBackSubstitute(temp,nRows,indx,col);
                setcolumn(result, j, col);
            }

        }
        else
        {
            int n = ((nRows)>=(nCols)?(nRows):(nCols));
            double temp[][] = identitymat(n);
            copymat(temp, M);



            invertmat(temp, temp);

            copymat(result, temp);
        }
    }

    public static double detDestructive(double M[][])
    {
        int n = M.length;
        { if (!(n == ((M).length==0 ? 0 : (M)[0].length))) HyperbolicUtils.set_exception("VecMath.java","detDestructive" , "n == ROWLENGTH(M)" ); };


        if (true)
        {

            switch (n)
            {
                case 0:
                    return 1.;
                case 1:
                    return M[0][0];
                case 2:
                    return vxv2(M[0], M[1]);
                case 3:
                    return vxvxv3(M[0], M[1], M[2]);
            }
        }
        else
        {
            if (((0 <=(n))&&((n)<= 3)))
                System.err.println("WARNING: doing slow detDestructive for testing, dim = "+n+"");
        }
        int indx[] = new int[n];
        try
        {
            double d = (double)luDecompose(M, n, indx);
            for (int j = 0; (j) < (n); ++j)
                d *= M[j][j];
            return d;
        }
        catch (Exception e)
        {
            return 0.;
        }
    }
    public static double det(double M[][])
    {
        int n = M.length;
        { if (!(n == ((M).length==0 ? 0 : (M)[0].length))) HyperbolicUtils.set_exception("VecMath.java","det" , "n == ROWLENGTH(M)" ); };


        if (true)
        {

            switch (n)
            {
                case 0:
                    return 1.;
                case 1:
                    return M[0][0];
                case 2:
                    return vxv2(M[0], M[1]);
                case 3:
                    return vxvxv3(M[0], M[1], M[2]);
            }
        }
        else
        {
            if (((0 <=(n))&&((n)<= 3)))
                System.err.println("WARNING: doing slow det for testing, dim = "+n+"");
        }

        return detDestructive(copymat(M));
    }


    public static void invmxv(double result[], double M[][], double b[])
    {
        int n = result.length;
        { if (!(n == M.length)) HyperbolicUtils.set_exception("VecMath.java","invmxv" , "n == M.length" ); };
        { if (!(n == ((M).length==0 ? 0 : (M)[0].length))) HyperbolicUtils.set_exception("VecMath.java","invmxv" , "n == ROWLENGTH(M)" ); };

        double temp[][] = copymat(M);
        int indx[] = new int[n];
        try
        {
                      luDecompose(temp, n, indx);
        }
        catch (Exception e)
        {
            fillvec(result, Double.POSITIVE_INFINITY);
            return;
        }
        copyvec(result, b);
        luBackSubstitute(temp, n, indx, result);

        if (normsqrd(result) < 1e20)
        {

            { if (!(result != b)) HyperbolicUtils.set_exception("VecMath.java","invmxv" , "result != b" ); };
            double bb[] = mxv(M, result);





            for (int i = 0; (i) < (n); ++i)
            {

                if (1. + ((bb[i] - b[i])*(bb[i] - b[i])) != 1.)
                {
                    {System.out.print("    ");System.out.println("\"WARNING!\"" + " = " + ("WARNING!"));};
                    {System.out.print("    ");System.out.println("bb[i]" + " = " + (bb[i]));};
                    {System.out.print("    ");System.out.println("b[i]" + " = " + (b[i]));};
                    {System.out.print("    ");System.out.println("((bb[i] - b[i])*(bb[i] - b[i]))" + " = " + (((bb[i] - b[i])*(bb[i] - b[i]))));};
                    {System.out.print("    ");System.out.println("1. + ((bb[i] - b[i])*(bb[i] - b[i]))" + " = " + (1. + ((bb[i] - b[i])*(bb[i] - b[i]))));};
                }
                { if (!((float)(1. + ((bb[i] - b[i])*(bb[i] - b[i]))) == (float)1.)) HyperbolicUtils.set_exception("VecMath.java","invmxv" , "(float)(1. + SQR(bb[i] - b[i])) == (float)1." ); };
            }
        }
    }

    public void invmxm(double result[][], double M[][], double B[][])
    {
        int n = result.length;
        { if (!(n == B.length)) HyperbolicUtils.set_exception("VecMath.java","invmxm" , "n == B.length" ); };
        { if (!(n == M.length)) HyperbolicUtils.set_exception("VecMath.java","invmxm" , "n == M.length" ); };
        { if (!(n == ((M).length==0 ? 0 : (M)[0].length))) HyperbolicUtils.set_exception("VecMath.java","invmxm" , "n == ROWLENGTH(M)" ); };
        int nCols = ((result).length==0 ? 0 : (result)[0].length);
        { if (!(nCols == ((B).length==0 ? 0 : (B)[0].length))) HyperbolicUtils.set_exception("VecMath.java","invmxm" , "nCols == ROWLENGTH(B)" ); };

        double temp[][] = copymat(M);
        int indx[] = new int[n];
        try
        {
                      luDecompose(temp, n, indx);
        }
        catch (Exception e)
        {
            fillmat(result, Double.POSITIVE_INFINITY);
            return;
        }
        double col[] = new double[n];
        for (int i = 0; (i) < (nCols); ++i)
        {
            getcolumn(col, B, i);
            luBackSubstitute(temp, n, indx, col);
            setcolumn(result, i, col);
        }
    }

    public void vxinvm(double result[], double v[], double M[][])
    {
        invmxv(result, transpose(M), v);
    }

    public void mxinvm(double result[][], double v[][], double M[][])
    {
        double transposeResult[][] = new double[((result).length==0 ? 0 : (result)[0].length)][result.length];
        invmxm(transposeResult, transpose(M), transpose(v));
        transpose(result, transposeResult);
    }

    public static void submat(double result[][], double M[][], int inds[])
    {
        for (int i = (inds.length)-1; (i) >= 0; --i)
        for (int j = (inds.length)-1; (j) >= 0; --j)
            result[i][j] = M[inds[i]][inds[j]];
    }

    public static void identitymat(double M[][])
    {
        for (int i = (M.length)-1; (i) >= 0; --i)
        for (int j = (M[i].length)-1; (j) >= 0; --j)
            M[i][j] = (i==j ? 1. : 0);
    }

    public static void makeRowRotMat(double M[][], int fromAxis, int toAxis, double radians)
    {
        identitymat(M);
        double s = Math.sin(radians);
        double c = Math.cos(radians);
        M[fromAxis][fromAxis] = c;
        M[fromAxis][toAxis] = s;
        M[toAxis][fromAxis] = -s;
        M[toAxis][toAxis] = c;
    }

    public static void makeRowTransMat(double M[][], double translate[])
    {
        identitymat(M);
        int nRows = M.length;
        for (int i = (translate.length)-1; (i) >= 0; --i)
            M[nRows-1][i] = translate[i];
    }
    public static void makeRowTransMatInv(double M[][], double translate[])
    {
        identitymat(M);
        int nRows = M.length;
        for (int i = (translate.length)-1; (i) >= 0; --i)
            M[nRows-1][i] = -translate[i];
    }




    public static void makeRowScaleMat(double M[][], double scale[])
    {
        int i, j, n = M.length, m = (n == 0 ? 0 : M[0].length);
        for (i = 0; (i) < (n); ++i)
        for (j = 0; (j) < (m); ++j)
            M[i][j] = (i==j ? i>=scale.length ? 1.
                                              : scale[i]
                            : 0);
    }

    public static double angleBetweenUnitVectors(double u0[], double u1[])
    {

        if (dot(u0, u1) > 0)
            return 2*Math.asin(.5*dist(u0,u1));
        else
        {
            double minusU1[] = sxv(-1.,u1);
            return Math.PI - 2*Math.asin(.5*dist(u0,minusU1));
        }
    }

    public static double[][] insertIdentityRowAndColumn(double M[][], int ind, double zero, double one)
    {
        double[][] result = new double[M.length+1][M.length+1];
        int i, j;
        for (i = 0; (i) < (M.length); ++i)
        for (j = 0; (j) < (M.length); ++j)
            result[i>=ind?i+1:i][j>=ind?j+1:j] = M[i][j];
        for (i = 0; (i) < (result.length); ++i)
            result[i][ind] = result[ind][i] = zero;
        result[ind][ind] = one;
        return result;
    }
    public static double[][] deleteRowAndColumn(double M[][], int ind)
    {
        double[][] result = new double[M.length-1][M.length-1];
        int i, j, n = result.length;
        for (i = 0; (i) < (n); ++i)
        for (j = 0; (j) < (n); ++j)
            result[i][j] = M[i>=ind?i+1:i][j>=ind?j+1:j];
        return result;
    }


    public static boolean equals(int v0[], int v1[])
    {
        if (v0 == null)
            return v1 == null ? true : false;
        else if (v1 == null)
            return false;
        int n = v0.length;
        if (n != v1.length)
            return false;
        for (int i = 0; (i) < (n); ++i)
            if (v0[i] != v1[i])
                return false;
        return true;
    }
    public static boolean equalsExactly(double v0[], double v1[])
    {
        if (v0 == null)
            return v1 == null ? true : false;
        else if (v1 == null)
            return false;
        int n = v0.length;
        if (n != v1.length)
            return false;
        for (int i = 0; (i) < (n); ++i)
            if (v0[i] != v1[i])
                return false;
        return true;
    }
    public static boolean equalsUsingEqualsSymbol(Object v0[], Object v1[])
    {
        if (v0 == null)
            return v1 == null ? true : false;
        else if (v1 == null)
            return false;
        int n = v0.length;
        if (n != v1.length)
            return false;
        for (int i = 0; (i) < (n); ++i)
            if (v0[i] != v1[i])
                return false;
        return true;
    }


    public static boolean equals(double M0[][], double M1[][], double eps)
    {
        if (M0 == null && M1 == null)
            return true;
        if (M0 == null || M1 == null || M0.length != M1.length)
            return false;
        int i, n = M0.length;
        for (i = 0; (i) < (n); ++i)
        {
            if (M0[i].length != M1[i].length)
                return false;
            int j, m = M0[i].length;
            for (j = 0; (j) < (m); ++j)
                if (M0[i][j] != M1[i][j]
                 && Math.abs(M0[i][j]-M1[i][j]) > eps)
                    return false;
        }
        return true;
    }


    public static void downshift(double to[][], double from[][])
    {
        double from0[] = from[0];
        int i, n = from.length;
        for (i = 0; (i) < (n-1); ++i)
            to[i] = from[i+1];
        to[i] = from0;
    }


    public static void upshift(double to[][], double from[][])
    {
        int i, n = from.length;
        double fromLast[] = from[n-1];
        for (i = n-1; i > 0; --i)
            to[i] = from[i-1];
        to[i] = fromLast;
    }

        public static int jacobi(double A[][], double d[], double V[][])
        {
            int n = d.length;

            double b[] = new double[n];
            double z[] = new double[n];
            identitymat(V);


            {
                int i;
                for (i = 0; (i) < (n); ++i)
                {
                    b[i] = d[i] = A[i][i];
                    z[i] = 0.;
                }
            }

            int nRot = 0;
            int iIter;
            for (iIter = 0; (iIter) < (50); ++iIter)
            {
                double sum = 0.;
                {
                    int i, j;
                    for (i = 0; (i) < (n-1); ++i)
                    for (j = i+1; j < n; ++j)
                        sum += ((A[i][j]) < 0 ? -(A[i][j]) : (A[i][j]));
                }

                if (sum == 0.)
                {



                    return nRot;
                }
                double tresh = (iIter < 3 ? .2*sum/(n*n) : 0.);
                int i, j;
                for (i = 0; (i) < (n-1); ++i)
                for (j = i+1; j < n; ++j)
                {
                    double g = 100.*((A[i][j]) < 0 ? -(A[i][j]) : (A[i][j]));

                    if (iIter > 3 && ((d[i]) < 0 ? -(d[i]) : (d[i]))+g == ((d[i]) < 0 ? -(d[i]) : (d[i]))
                                  && ((d[j]) < 0 ? -(d[j]) : (d[j]))+g == ((d[j]) < 0 ? -(d[j]) : (d[j])))
                        A[i][j] = 0.;
                    else if (((A[i][j]) < 0 ? -(A[i][j]) : (A[i][j])) > tresh)
                    {
                        double h = d[j] - d[i];
                        double t;
                        if (((h) < 0 ? -(h) : (h))+g == ((h) < 0 ? -(h) : (h)))
                        {
                            t = A[i][j]/h;
                        }
                        else
                        {
                            double theta = .5*h/(A[i][j]);
                            t = 1./(((theta) < 0 ? -(theta) : (theta)) + Math.sqrt(1. + theta*theta));
                            if (theta < 0.)
                                t = -t;
                        }
                        double c = 1./Math.sqrt(1. + t*t);
                        double s = t*c;
                        double tau = s/(1. + c);
                        h = t * A[i][j];
                        z[i] -= h;
                        z[j] += h;
                        d[i] -= h;
                        d[j] += h;
                        A[i][j] = 0.;
                        int k;

                        for (k = 0; (k) < (i); ++k)
                            { g = A[k][i]; h = A[k][j]; A[k][i] = g-s*(h+g*tau); A[k][j] = h+s*(g-h*tau); };
                        for (k = i+1; k < j; ++k)
                            { g = A[i][k]; h = A[k][j]; A[i][k] = g-s*(h+g*tau); A[k][j] = h+s*(g-h*tau); };
                        for (k = j+1; k < n; ++k)
                            { g = A[i][k]; h = A[j][k]; A[i][k] = g-s*(h+g*tau); A[j][k] = h+s*(g-h*tau); };
                        for (k = 0; (k) < (n); ++k)
                            { g = V[k][i]; h = V[k][j]; V[k][i] = g-s*(h+g*tau); V[k][j] = h+s*(g-h*tau); };
                        nRot++;
                    }
                }

                for (i = 0; (i) < (n); ++i)
                {
                    b[i] += z[i];
                    d[i] = b[i];
                    z[i] = 0.;
                }
            }



            System.err.println("Too many iterations in routine jacobi");
            return -1;
        }


    public static int leftSqrtOfSymmetricMatrix(double result[][], boolean whichColumnsOfResultAreImaginary[], double C[][])
    {
        int n = C.length;
        { if (!(n == result.length)) HyperbolicUtils.set_exception("VecMath.java","leftSqrtOfSymmetricMatrix" , "n == result.length" ); };
        double eigenValues[] = new double[n];
        double eigenVectors[][] = new double[n][n];
        double temp[][] = copymat(C);
        int nRot = jacobi(temp, eigenValues, eigenVectors);
        if (nRot < 0)
        {


            if (result != null)
                identitymat(result);
            return -2;
        }

        if (true)
        {


            double diagEigenValues[][] = identitymat(n);
            int i;
            for (i = 0; (i) < (n); ++i)
                diagEigenValues[i][i] = eigenValues[i];
            double A[][] = mxm(C, eigenVectors);
            double B[][] = mxm(eigenVectors, diagEigenValues);





            int j;
            for (i = 0; (i) < (n); ++i)
            for (j = 0; (j) < (n); ++j)
                { if (!(1. + ((A[i][j]-B[i][j])*(A[i][j]-B[i][j])) == 1.)) HyperbolicUtils.set_exception("VecMath.java","leftSqrtOfSymmetricMatrix" , "1. + SQR(A[i][j]-B[i][j]) == 1." ); };
        }

        int positiveDefiniteStatus = 1;
        int i;
        for (i = 0; (i) < (n); ++i)
        {
            double eigenValue = eigenValues[i];
            if (1.+((eigenValue)*(eigenValue)) == 1.)
                positiveDefiniteStatus = ((positiveDefiniteStatus)<=(0)?(positiveDefiniteStatus):(0));
            else if (eigenValue < 0.)
                positiveDefiniteStatus = -1;

            double sqrtEigenValue;
            if (whichColumnsOfResultAreImaginary != null)
            {
                if ((whichColumnsOfResultAreImaginary[i] = (eigenValue < 0.)) == true)
                    sqrtEigenValue = Math.sqrt(-eigenValue);
                else
                    sqrtEigenValue = Math.sqrt(eigenValue);
            }
            else
                sqrtEigenValue = Math.sqrt(((eigenValue)>=(0.)?(eigenValue):(0.)));


            int j;
            for (j = 0; (j) < (n); ++j)
                eigenVectors[j][i] *= sqrtEigenValue;
        }

        copymat(result, eigenVectors);
        return positiveDefiniteStatus;
    }

    public static int positiveSemiDefinitizeWithUnitDiagonal(double result[][], boolean whichColumnsAreImaginary[], double C[][])
    {
        double B[][] = new double[result.length][result.length];
        int positiveDefiniteStatus = leftSqrtOfSymmetricMatrix(B, whichColumnsAreImaginary, C);
        if (positiveDefiniteStatus == -2)
        {

            identitymat(result);
        }
        else
        {
            if (whichColumnsAreImaginary == null)
                correlate(result, B);
        }
        return positiveDefiniteStatus;
    }

    public static int positiveDefiniteness(double M[][])
    {
        return positiveSemiDefinitizeWithUnitDiagonal(null, null, M);
    }


    public static void correlate(double result[][], double V[][])
    {
        double normalizedV[][] = new double[V.length][V.length==0?0:V[0].length];
        normalizeRows(normalizedV, V);
        mxm(result, normalizedV, transpose(normalizedV));
    }

    public static void normalizeRows(double result[][], double M[][])
    {
        int i, n = M.length;
        for (i = 0; (i) < (n); ++i)
            normalize(result[i], M[i]);
    }


    public static void normalize(double result[], double v[])
    {
        double len = norm(v);
        vxs(result, v, len==0. ? 1. : 1./len);
    }

    public static String toString(Object v)
    {
        if (v == null)
            return "(null)";
        int n = java.lang.reflect.Array.getLength(v);
        StringBuffer sb = new StringBuffer();
        sb.append("<");
        for (int i = 0; (i) < (n); ++i)
        {
            sb.append(java.lang.reflect.Array.get(v, i));
            if (i < n-1)
                sb.append(",");
        }
        sb.append(">");
        return sb.toString();
    }

    public static String toString(double m[][])
    {
        return Arrays.toStringNonCompact(m, "", "    ");
    }

    public static String toString(int m[][])
    {
        return Arrays.toStringNonCompact(m, "", "    ");
    }


    public static String toString(int m[][][], String bigSep, String littleSep)
    {
        if (m == null)
            return "(null)";
        String result = "";
        int iLayer, nLayers = m.length;
        for (iLayer = 0; (iLayer) < (nLayers); ++iLayer)
        {
            int iRow, nRows = m[iLayer].length;
            for (iRow = 0; (iRow) < (nRows); ++iRow)
            {
                result += "    [";
                int iCol, nCols = m[iLayer][iRow].length;
                for (iCol = 0; (iCol) < (nCols); ++iCol)
                {
                    result += m[iLayer][iRow][iCol];
                    if (iCol+1 < nCols)
                        result += " ";
                }
                result += "]";
                if (iRow+1 < nRows)
                    result += littleSep;
            }
            if (iLayer+1 < nLayers)
                result += bigSep;
        }
        return result;
    }
    public static String toString(int m[][][])
    {
        return toString(m, "\n\n", "\n");

    }







        public static double[] vpv(double v0[], double v1[])
        {
            double result[] = new double[v0.length];
            vpv(result, v0, v1);
            return result;
        }
        public static int[] vpv(int v0[], int v1[])
        {
            int result[] = new int[v0.length];
            vpv(result, v0, v1);
            return result;
        }
        public static double[] sxv(double s, double v[])
        {
            double result[] = new double[v.length];
            sxv(result, s, v);
            return result;
        }
        public static double[] vmv(double v0[], double v1[])
        {
            double result[] = new double[v0.length];
            vmv(result, v0, v1);
            return result;
        }
        public static double[][] mmv(double M[][], double v[])
        {
            double result[][] = new double[M.length][v.length];
            mmv(result, M, v);
            return result;
        }
        public static double[] lerp(double p0[], double p1[], double t)
        {
            double result[] = new double[p0.length];
            lerp(result, p0, p1, t);
            return result;
        }
        public static double[] bary(double p0[],
                                    double p1[], double t1,
                                    double p2[], double t2)
        {
            double result[] = new double[p0.length];
            bary(result, p0, p1, t1, p2, t2);
            return result;
        }
        public static double[] zerovec(int dim)
        {
            double result[] = new double[dim];
            zerovec(result);
            return result;
        }
        public static double[] fillvec(int dim, double s)
        {
            double result[] = new double[dim];
            fillvec(result,s);
            return result;
        }
        public static int[] fillvec(int dim, int s)
        {
            int result[] = new int[dim];
            fillvec(result,s);
            return result;
        }
        public static boolean[] fillvec(int dim, boolean s)
        {
            boolean result[] = new boolean[dim];
            fillvec(result,s);
            return result;
        }
        public static double[][] fillmat(int dim0, int dim1, double s)
        {
            double result[][] = new double[dim0][dim1];
            fillmat(result,s);
            return result;
        }
        public static double[] sum(double array[][])
        {
            if (array.length == 0)
                return null;
            double result[] = new double[array[0].length];
            sum(result, array);
            return result;
        }
        public static double[] average(double array[][])
        {
            if (array.length == 0)
                return null;
            double result[] = new double[array[0].length];
            average(result, array);
            return result;
        }
        public static double[] averageIndexed(int inds[][], double array[][])
        {
            if (array.length == 0)
                return null;
            double result[] = new double[array[0].length];
            averageIndexed(result, inds, array);
            return result;
        }
        public static double[ ][] bboxIndexed(double array[][],
                                                  Object inds)
        {
            if (array.length == 0)
                return null;
            double result[][] = new double[2][array[0].length];
            bboxIndexed(result, array, inds);
            return result;
        }
        public static double[ ][] bboxUniform(double array[][])
        {
            if (array.length == 0)
                return null;
            double result[][] = new double[2][array[0].length];
            bboxUniform(result, array);
            return result;
        }
        public static double[ ][] bbox(double array[][])
        {
            if (array.length == 0)
                return null;
            double result[][] = new double[2][array[0].length];
            bbox(result, array);
            return result;
        }
        public static double[ ][] bboxIntersect(double bbox0[ ][],
                                                    double bbox1[ ][])
        {
            double result[][] = new double[2][bbox0[0].length];
            bboxIntersect(result, bbox0, bbox1);
            return result;
        }
        public static double[ ][] bboxUnion(double bbox0[ ][],
                                                double bbox1[ ][])
        {
            double result[][] = new double[2][bbox0[0].length];
            bboxUnion(result, bbox0, bbox1);
            return result;
        }

        public static double[][] submat(double M[][], int inds[])
        {
            double result[][] = new double[inds.length][inds.length];
            submat(result, M, inds);
            return result;
        }
        public static double[][] identitymat(int n)
        {
            double result[][] = new double[n][n];
            identitymat(result);
            return result;
        }
        public static double[][] identitymat(int n, int m)
        {
            double result[][] = new double[n][m];
            identitymat(result);
            return result;
        }
        public static double[][] transpose(double M[][])
        {
            double result[][] = new double[M.length==0 ? 0 : M[0].length][M.length];
            transpose(result, M);
            return result;
        }
        public static double[][] copymat(double M[][])
        {
            double result[][] = new double[M.length][M.length==0 ? 0 : M[0].length];
            copymat(result, M);
            return result;
        }
        public static boolean[] copyvec(boolean v[])
        {
            boolean result[] = new boolean[v.length];
            copyvec(result, v);
            return result;
        }
        public static int[] copyvec(int v[])
        {
            int result[] = new int[v.length];
            copyvec(result, v);
            return result;
        }
        public static double[] copyvec(double v[])
        {
            double result[] = new double[v.length];
            copyvec(result, v);
            return result;
        }
        public static double[] getcolumn(double M[][], int iCol)
        {
            double result[] = new double[M.length];
            getcolumn(result, M, iCol);
            return result;
        }
        public static double[][] mxm(double m0[ ][ ],
                                     double m1[ ][ ])
        {

            int n = m0.length;
            int dotLength = m1.length;
            int m = (dotLength == 0 ? 0 : m1[0].length);

            {

                if (dotLength != ((m0).length==0 ? 0 : (m0)[0].length))
                {
                    { if (!(dotLength == 4 && m0[0].length == 3)) HyperbolicUtils.set_exception("VecMath.java","mxm" , "dotLength == 4 && m0[0].length == 3" ); };
                    n = 4;
                    m = 3;
                }
            }

            double[][] result = new double[n][m];
            mxm(result, m0, m1);
            return result;

        }

        public static int[][] mxm(int m0[ ][ ],
                                  int m1[ ][ ])
        {
            int n = m0.length;
            int dotLength = m1.length;
            int m = (dotLength == 0 ? 0 : m1[0].length);
            int[][] result = new int[n][m];
            mxm(result, m0, m1);
            return result;
        }
        public static double[] sxvpsxv(double s0, double v0[],
                                       double s1, double v1[])
        {
            double[] result = new double[v0.length];
            sxvpsxv(result, s0, v0, s1, v1);
            return result;
        }
        public static double[][] mxs(double M[][], double s)
        {
            double result[][] = new double[M.length][M.length==0 ? 0 : M[0].length];
            mxs(result, M, s);
            return result;
        }
        public static double[] vxs(double v[], double s)
        {
            double result[] = new double[v.length];
            vxs(result, v, s);
            return result;
        }
        public static double[] vxm(double v[], double M[][])
        {
            double result[] = new double[M.length == 0 ? 0 : M[0].length];
            vxm(result, v, M);
            return result;
        }
        public static double[] mxv(double M[][], double v[])
        {
            double result[] = new double[M.length];
            mxv(result, M, v);
            return result;
        }
        public static double[] vpsxv(double v0[], double s1, double v1[])
        {
            double result[] = new double[v0.length];
            vpsxv(result, v0, s1, v1);
            return result;
        }
        public static int[] vpsxv(int v0[], int s1, int v1[])
        {
            int result[] = new int[v0.length];
            vpsxv(result, v0, s1, v1);
            return result;
        }
        public static double[][] invertmat(double M[][])
        {
            double result[][] = new double[M.length][M.length==0 ? 0 : M[0].length];
            invertmat(result, M);
            return result;
        }
        public static double[] invmxv(double M[][], double v[])
        {
            double result[] = new double[M.length];
            invmxv(result, M, v);
            return result;
        }

        public static double[] vxv3(final double v1[ ], final double v2[ ])
        {
            double result[] = new double[3];
            vxv3(result, v1, v2);
            return result;
        }

        public static double[] crossprod(double vectors[][])
        {
            double result[] = new double[vectors.length+1];
            crossprod(result, vectors);
            return result;
        }

        public static double[] random(int dim)
        {
            double result[] = new double[dim];
            random(result);
            return result;
        }
        public static double[] random(int dim, java.util.Random generator)
        {
            double result[] = new double[dim];
            random(result, generator);
            return result;
        }


        public static double[][] makeRowTransMat(double translate[])
        {
            double result[][] = new double[translate.length+1][translate.length];
            makeRowTransMat(result, translate);
            return result;
        }

        public static double[][] makeRowTransMatInv(double translate[])
        {
            double result[][] = new double[translate.length+1][translate.length];
            makeRowTransMatInv(result, translate);
            return result;
        }

        public static double[][] makeRowTransMat(double x, double y, double z)
        {
            return makeRowTransMat(new double[] {x, y, z});
        }

        public static double[][] makeRowTransMatInv(double x, double y, double z)
        {
            return makeRowTransMatInv(new double[] {x, y, z});
        }

        public static double[][] makeRowScaleMat(double scale[])
        {
            double result[][] = new double[scale.length+1][scale.length];
            makeRowScaleMat(result, scale);
            return result;
        }
        public static double[][] makeRowScaleMat(double scale[],
                        double tiePointIn[],
                        double tiePointOut[])
        {
            double M[][] = makeRowScaleMat(scale);
            if (tiePointIn != null)
                M = mxm(makeRowTransMatInv(tiePointIn), M);
            if (tiePointOut != null)
                M = mxm(M, makeRowTransMat(tiePointOut));
            return M;
        }
        public static double[][] makeRowScaleMat(double scale[],
                        double fixedPoint[])
        {
            return makeRowScaleMat(scale, fixedPoint, fixedPoint);
        }
        public static double[][] makeRowScaleMat(int n, double scale,
                        double tiePointIn[],
                        double tiePointOut[])
        {
            return makeRowScaleMat(fillvec(n,scale),tiePointIn,tiePointOut);
        }
        public static double[][] makeRowScaleMat(int n, double scale,
                                                 double fixedPoint[])
        {
            return makeRowScaleMat(n,scale,fixedPoint,fixedPoint);
        }
        public static double[][] makeRowScaleMat(int n, double scale)
        {
            return makeRowScaleMat(n,scale,null,null);
        }
        public static double[][] makeRowScaleMat(double x, double y, double z)
        {
            return makeRowScaleMat(new double[]{x,y,z},null,null);
        }




        public static double[][] makeRowRotMat(int n, int fromAxis, int toAxis, double radians)
        {
            double result[][] = new double[n][n];
            makeRowRotMat(result, fromAxis, toAxis, radians);
            return result;
        }

        public static double[][] makeRowRotMat(int n, int fromAxis, int toAxis, double radians,
                double tiePointIn[],
                double tiePointOut[])
        {
            double M[][] = makeRowRotMat(n, fromAxis, toAxis, radians);
            if (tiePointIn != null)
                M = mxm(makeRowTransMatInv(tiePointIn), M);
            if (tiePointOut != null)
                M = mxm(M, makeRowTransMat(tiePointOut));
            return M;
        }

        public static double[][] makeRowRotMat(int n, int fromAxis, int toAxis, double radians, double fixedPoint[])
        {
            return makeRowRotMat(n, fromAxis, toAxis, radians, fixedPoint, fixedPoint);
        }

        public static double[][] makeRowRotMat(double radians, double fixedAxes[ ][])
        {
            int nFixedAxes = fixedAxes.length;
            int n = (nFixedAxes == 0 ? 2 : fixedAxes[0].length);
            { if (!(nFixedAxes == n-2)) HyperbolicUtils.set_exception("VecMath.java","makeRowRotMat" , "nFixedAxes == n-2" ); };

            { if (!(n == 3)) HyperbolicUtils.set_exception("VecMath.java","makeRowRotMat" , "n == 3" ); };

            double fixedAxis[] = fixedAxes[0];

            double rotPlusZToFixedAxis[][];
            {
                rotPlusZToFixedAxis = new double[3][3];
                VecMath.copyvec(rotPlusZToFixedAxis[2], fixedAxis);
                VecMath.normalize(rotPlusZToFixedAxis[2], rotPlusZToFixedAxis[2]);
                int mini = VecMath.mini(fixedAxis);
                VecMath.zerovec(rotPlusZToFixedAxis[1]);
                rotPlusZToFixedAxis[1][mini] = 1.;

                VecMath.vxv3(rotPlusZToFixedAxis[0],
                             rotPlusZToFixedAxis[1],
                             rotPlusZToFixedAxis[2]);
                VecMath.normalize(rotPlusZToFixedAxis[0], rotPlusZToFixedAxis[0]);
                VecMath.vxv3(rotPlusZToFixedAxis[1],
                             rotPlusZToFixedAxis[2],
                             rotPlusZToFixedAxis[0]);
            }

            return mxmxm(VecMath.transpose(rotPlusZToFixedAxis),
                         makeRowRotMat(3, 0, 1, radians),
                         rotPlusZToFixedAxis);
        }

        public static double[ ][ ] makeRowTiePointMat(
                               double inTiePoints[][],
                               double outTiePoints[][])
        {
            { if (!(inTiePoints.length > 0)) HyperbolicUtils.set_exception("VecMath.java","makeRowTiePointMat" , "inTiePoints.length > 0" ); };

            int dim = inTiePoints[0].length;
            { if (!(inTiePoints.length == dim+1)) HyperbolicUtils.set_exception("VecMath.java","makeRowTiePointMat" , "inTiePoints.length == dim+1" ); };
            { if (!(outTiePoints[0].length == dim)) HyperbolicUtils.set_exception("VecMath.java","makeRowTiePointMat" , "outTiePoints[0].length == dim" ); };
            { if (!(outTiePoints.length == dim+1)) HyperbolicUtils.set_exception("VecMath.java","makeRowTiePointMat" , "outTiePoints.length == dim+1" ); };



            double inPadded[][] = VecMath.fillmat(dim+1,dim+1, 1.);
            VecMath.copymat(inPadded, inTiePoints);
            double Mpadded[][] = VecMath.mxm(VecMath.invertmat(inPadded),
                                             outTiePoints);
            double M[][] = new double[dim+1][dim];
            VecMath.copymat(M, Mpadded);
            return M;
        }



        public static double[] flatten(double M[][])
        {
            return (double[])Arrays.flatten(M, 0, 2);
        }

        public static double[][] mxmxm(double A[][], double B[][], double C[][])
        {
            return mxm(mxm(A,B),C);
        }
        public static double[][] mxmxmxm(double A[][], double B[][], double C[][], double D[][])
        {
            return mxm(mxmxm(A,B,C),D);
        }
        public static double[] vxmxm(double a[], double B[][], double C[][])
        {

            return vxm(vxm(a,B),C);
        }


        public static int permutationSignDestructive(int perm[])
        {
            int sign = 1;
            int permi;
            int n = perm.length;
            for (int i = 0; (i) < (n); ++i)
            {
                permi = perm[i];
                while (permi != i)
                {

                    {


                        int temp = perm[permi];
                        perm[permi] = permi;
                        permi = temp;
                    }
                    sign = -sign;
                }
            }
            return sign;
        }

        public static boolean isIdentityPerm(int perm[])
        {
            for (int i = (perm.length)-1; (i) >= 0; --i)
                if (perm[i] != i)
                    return false;
            return true;
        }
        public static void identityperm(int perm[])
        {
            for (int i = (perm.length)-1; (i) >= 0; --i)
                perm[i] = i;
        }
        public static int[] identityperm(int n)
        {
            int result[] = new int[n];
            identityperm(result);
            return result;
        }
        public static void randomperm(int perm[], java.util.Random generator)
        {
            identityperm(perm);
            Arrays.shuffle(perm, generator);
        }
        public static int[] randomperm(int n, java.util.Random generator)
        {
            int result[] = new int[n];
            randomperm(result, generator);
            return result;
        }

        public static void invertperm(int result[], int perm[])
        {
            { if (!(result != perm)) HyperbolicUtils.set_exception("VecMath.java","invertperm" , "result != perm" ); };
            int n = ((result.length)<=(perm.length)?(result.length):(perm.length));
            for (int i = (n)-1; (i) >= 0; --i)
                result[perm[i]] = i;
        }
        public static int[] invertperm(int perm[])
        {
            int result[] = new int[perm.length];
            invertperm(result, perm);
            return result;
        }






        public static void composeluts(int result[],
                                       int lut0[],
                                       int lut1[])
        {
            int n = lut1.length;
            { if (!(n == result.length)) HyperbolicUtils.set_exception("VecMath.java","composeluts" , "n == result.length" ); };
            for (int i = 0; (i) < (n); ++i)
                result[i] = lut0[lut1[i]];
        }





        public static int[] composeluts(int lut0[],
                                        int lut1[])
        {
            int result[] = new int[lut1.length];
            composeluts(result, lut0, lut1);
            return result;
        }

    static void main(String args[])
    {
        if (args.length < 1)
        {
            System.err.println("Usage: VecMath <maxdim> [<expr>]");
            System.exit(1);
        }

        int mindim = 2;
        int maxdim = Integer.parseInt(args[0]);

        for (int dim = mindim; dim <= maxdim; ++dim)
        {
            {
                double M[][] = new double[dim][dim];
                double a[] = new double[dim];
                double b[] = new double[dim];
                double x[] = new double[dim];
                for (int i = 0; (i) < (dim); ++i)
                    random(M[i]);
                random(a);


                mxv(b,M,a);
                invmxv(x,M,b);

                System.out.println("    " + "M" + " =\n" + VecMath.toString(M));
                System.out.println("    " + "a" + " = " + VecMath.toString(a));
                System.out.println("    " + "b" + " = " + VecMath.toString(b));
                System.out.println("    " + "x" + " = " + VecMath.toString(x));
            }

            if (dim >= 1)
            {
                double vectors[][] = new double[dim-1][dim];
                identitymat(vectors);
                System.out.println("    " + "vectors" + " =\n" + VecMath.toString(vectors));
                System.out.println("    " + "crossprod(vectors)" + " = " + VecMath.toString(crossprod(vectors)));
            }
        }

    }

}
