import javax.swing.JOptionPane;

public class HyperbolicUtils
{
	public static void message(String str)
		{
		JOptionPane.showMessageDialog(null,str);
		}
	public static void set_exception(String fichier, String procedure, String libelle)
		{
		String str=libelle+" - fichier : "+fichier+" procedure : "+procedure;
		JOptionPane.showMessageDialog(null,str);
		throw new Error(str);
		}

    public static double
    h2eNorm(double hNorm)
    {
        if (Double.isInfinite(hNorm))
            return 1.;
        return MyMath.tanh(.5*hNorm);
    }



    public static double
    e2hNorm(double eNorm)
    {
        return 2*MyMath.atanh(eNorm);
    }

    public static boolean
    unHBary2(int nVerts,
             double result[ ],
             double p[ ],
             double verts[ ][ ],
             int verbose)
    {
        final int dim = 2;
        { if (!(nVerts >= 3)) set_exception("HyperbolicUtils.java","unHBary2" , "nVerts >= 3"); };


        Isometry2 translatePToOrigin = Isometry2.pureTranslation(-p[0],-p[1]);
        double v[][] = new double[nVerts][dim];
        int i;
        for (i = 0; (i) < (nVerts); ++i)
        {
            Complex temp = translatePToOrigin.apply(new Complex(verts[i]));
            v[i][0] = temp.x;
            v[i][1] = temp.y;
        }

        for (i = 0; (i) < (3); ++i)
        {
            result[i] = calcSmallestPointOnLine2(null,
                                                 v[(i+1)%3],
                                                 v[(i+2)%3]);
            if (VecMath.vxv2(v[(i+1)%3],v[(i+2)%3]) < 0.)
                result[i] *= -1.;
        }

        for (; i < nVerts; ++i)
            result[i] = 0;

        return true;
    }


    public static boolean
    unHBary2(int nVerts,
             double result[],
             Complex p,
             Complex verts[],
             int verbose)
    {
        double _p[] = {p.x,p.y};
        double _verts[][] = new double[nVerts][2];
        int i;
        for (i = 0; (i) < (nVerts); ++i)
        {
            _verts[i][0] = verts[i].x;
            _verts[i][1] = verts[i].y;
        }
        return unHBary2(nVerts,
                        result,
                        _p,
                        _verts,
                        verbose);
    }

    public static boolean
    hBary2(int nPoints,
           double result[ ],
           double b[ ],
           Complex p[ ],
           int verbose)
    {
        double _p[][] = new double[nPoints][2];
        for (int i = 0; (i) < (nPoints); ++i)
        {
            _p[i][0] = p[i].x;
            _p[i][1] = p[i].y;
        }
        return hBary2(nPoints, result, b, _p, verbose);
    }




    public static boolean
    hBary2FromDistances(int nPoints,
                        double result[ ],
                        double dists[ ],
                        double p[ ][ ])
    {

        { if (!(nPoints >= 3)) set_exception("HyperbolicUtils.java","hBary2FromDistances","nPoints >= 3" ); };

        double p0[ ] = p[0];
        double p1[ ] = p[1];

        double A = dists[1];
        double B = dists[0];
        double C = hdist2_poincare(p0, p1);
        double alpha = solveTriangleAlphaFromABC(A, B, C);
        double beta = solveTriangleAlphaFromABC(B, C, A);

        if (pointIsToRightSideOfLinePoincare(p[2],
                                             p[0], p[1],
                                             0.))
        {
            { if (!(false)) set_exception("HyperbolicUtils.java","hBary2FromDistances" , "false"); };
            alpha *= -1.;
            beta *= -1.;
        }

        double foo0[];
        {
            Isometry2 take_p0_to_origin = Isometry2.pureTranslation(-p0[0], -p0[1]);
            Isometry2 take_origin_to_p0 = Isometry2.pureTranslation(p0[0], p0[1]);
            double q1[] = take_p0_to_origin.apply(p1);
            double dir = Math.atan2(q1[1], q1[0]) + alpha;
            foo0 = new double[] {.5*Math.cos(dir), .5*Math.sin(dir)};
            take_origin_to_p0.apply(foo0, foo0);
        }
        double foo1[];
        {
            Isometry2 take_p1_to_origin = Isometry2.pureTranslation(-p1[0], -p1[1]);
            Isometry2 take_origin_to_p1 = Isometry2.pureTranslation(p1[0], p1[1]);
            double q0[] = take_p1_to_origin.apply(p0);
            double dir = Math.atan2(q0[1], q0[0]) - beta;
            foo1 = new double[] {.5*Math.cos(dir), .5*Math.sin(dir)};
            take_origin_to_p1.apply(foo1, foo1);
        }

        return calcLineLineIntersection2(result,
                                         p0, foo0,
                                         p1, foo1);
    }


    public static boolean
    hBary2(int nPoints,
           double result[ ],
           double b[ ],
           double p[ ][ ],
           int verbose)
    {



        { if (!(nPoints >= 3)) set_exception("HyperbolicUtils.java","hBary2","nPoints >= 3"); };
        if (nPoints > 3)
        {
            return hBary2(3,
                          result,
                          b,
                          p,
                          verbose);
        }



        final int dim = 2;

        if (verbose >= 2)
        {
            System.out.println("in hBary");
            int i, j;
            System.out.print("    b =");
            for (i = 0; (i) < (nPoints); ++i)
                System.out.print(" " + b[i]);
            System.out.println();
            System.out.print("    p = ");
            for (i = 0; (i) < (nPoints); ++i)
            {
                for (j = 0; (j) < (dim); ++j)
                {
                    System.out.print("" + p[i][j] + (j==dim-1 ? " " : ","));
                }
            }
            System.out.println();
        }

        boolean doSimpleWrongWay = false;
        if (doSimpleWrongWay)
        {



            VecMath.zerovec(result);

            double sum = 0.;
            int i;
            for (i = 0; (i) < (nPoints); ++i)
            {
                sum += b[i];
                if (sum != 0.)
                {

                    Isometry2 temp = Isometry2.pureTranslation(result[0],result[1]);
                    Isometry2 inverseTemp = Isometry2.pureTranslation(-result[0],-result[1]);
                    double p_i_[] = inverseTemp.apply(p[i]);
                    VecMath.vxs(result, p_i_, b[i]/sum);
                    temp.apply(result,result);
                }
                else
                {


                    { if (!(b[i] == 0.)) set_exception("HyperbolicUtils.java","hBary2" ,"b[i] == 0." ); };
                }
            }
            return true;
        }


        double q[][] = new double[nPoints][dim];

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

            double otherPoints[][] = new double[nPoints-1][];
            double otherBs[] = new double[nPoints-1];
            int j;
            for (j = 0; (j) < (nPoints-1); ++j)
            {
                otherPoints[j] = p[(i+1+j)%nPoints];
                otherBs[j] = b[(i+1+j)%nPoints];
            }
            _calcPointOnWeightedSector2(nPoints-1, q[i], p[i], otherPoints, otherBs, verbose);
        }






        double bestSum = -1;
        int bestI = -1, bestJ = -1;
        for (i = 0; (i) < (nPoints); ++i)
        {
            int j;
            for (j = 0; (j) < (i); ++j)
                {
                    double sum = ((b[i] + b[j] - 1.) < 0 ? -(b[i] + b[j] - 1.) : (b[i] + b[j] - 1.));
                    if (sum > bestSum)
                    {
                        bestI = i;
                        bestJ = j;
                        bestSum = sum;
                    }
                }
        }
        { if (!(bestI != -1 && bestJ != -1 && bestI != bestJ)) set_exception("HyperbolicUtils.java","hBary2" , "bestI != -1 && bestJ != -1 && bestI != bestJ" ); };
        boolean returnval = calcLineLineIntersection2(result, p[bestI], q[bestI], p[bestJ], q[bestJ]);


        boolean hBaryDebug = false;
        if (hBaryDebug)
        {

            int j;
            for (i = 0; (i) < (nPoints); ++i)
            for (j = 0; (j) < (nPoints); ++j)
                if (i != j)
                {
                    double intersection[] = new double[dim];
                    System.out.print(
                           (i==bestI && j==bestJ ? "(best)" : "      ")
                           + " Intersection "+i+" "+j);
                    if (calcLineLineIntersection2(intersection, p[i], q[i], p[j], q[j]))
                    {
                        System.out.print(" =");
                        int k;
                        for (k = 0; (k) < (intersection.length); ++k)
                            System.out.print(" "+intersection[k]);
                        System.out.println();
                    }
                    else
                        System.out.println(" bad");
                }
        }
        return returnval;
    }






    public static boolean
    calcLineLineIntersection2(double result[ ],
                              double p0[ ], double p1[ ],
                              double q0[ ], double q1[ ])
    {
        final int N = 2;


        double q[] = new double[N];
        calcSmallestPointOnLine2(q,q0,q1);
        Isometry2 F = Isometry2.pureTranslation(q[0],q[1]);
        Isometry2 inverseF = F.inverse();

        double p0_[] = inverseF.apply(p0);
        double p1_[] = inverseF.apply(p1);
        double q0_[] = inverseF.apply(q0);
        double q1_[] = inverseF.apply(q1);


        p0_ = _calcEndOfGeodesic2(p1_,p0_);
        p1_ = _calcEndOfGeodesic2(p0_,p1_);
        q0_ = _calcEndOfGeodesic2(q1_,q0_);





        double x0 = VecMath.dot(q0_,p0_);
        double x1 = VecMath.dot(q0_,p1_);
        double temp = 2*x0*x1 - VecMath.dot(p0_,p1_) + 1;
        if (temp*temp-(x0+x1)*(x0+x1) < 0.)
            return false;
        double x = (x0+x1)/(temp+Math.sqrt(temp*temp-(x0+x1)*(x0+x1)));
        VecMath.vxs(result, q0_, x);


        F.apply(result, result);
        return true;
    }

    public static double
    calcSmallestPointOnLine2(double result[ ], double p0[ ], double p1[ ])
    {
        final int N = 2;


        double q0[] = _calcEndOfGeodesic2(p1, p0);
        double q1[] = _calcEndOfGeodesic2(p0, p1);


        double q[] = new double[N];
        VecMath.sxvpsxv(q, .5, q0, .5, q1);

        double scale = 1. / (1 + .5 * Math.sqrt(VecMath.distsqrd(q0,q1)));

        if (result != null)
            VecMath.vxs(result, q, scale);
        return e2hNorm(scale * Math.sqrt(VecMath.dot(q,q)));
    }

    public static double
    calcClosestPointOnLine2(double result[ ], double p[ ], double l0[ ], double l1[ ])
    {

        Isometry2 F = Isometry2.pureTranslation(p[0], p[1]);
        Isometry2 inverseF = F.inverse();
        double _l0[] = inverseF.apply(l0);
        double _l1[] = inverseF.apply(l1);
        double dist = calcSmallestPointOnLine2(result, _l0, _l1);
        F.apply(result, result);
        return dist;
    }


    public static double
    calcClosestPointOnLine2(Complex result, Complex p, Complex l0, Complex l1)
    {
        double _result[] = new double[2];
        double dist = calcClosestPointOnLine2(_result,
                                              new double[] {p.x,p.y},
                                              new double[] {l0.x,l0.y},
                                              new double[] {l1.x,l1.y});
        result.set(_result[0], _result[1]);
        return dist;
    }





    private static double[]
    _calcEndOfGeodesic2(double p0[ ], double p1[ ])
    {
        final int N = 2;


        Isometry2 F = Isometry2.pureTranslation(p1[0],p1[1]);
        Isometry2 inverseF = F.inverse();
        double p0_[] = inverseF.apply(p0);
        double len = Math.sqrt(VecMath.dot(p0_,p0_));
        { if (!(len != 0.)) set_exception("HyperbolicUtils.java","_calcEndOfGeodesic2" , "len != 0." ); };
        double result[] = new double[N];
        VecMath.vxs(result, p0_, -1./len);
        F.apply(result, result);
        return result;
    }

    private static void
    _calcPointOnWeightedSector2(int nPoints,
                                double result[ ],
                                double p0[ ],
                                double p[ ][ ],
                                double b[ ],
                                int verbose)
    {
        final int dim = 2;

        double d00 = VecMath.dot(p0,p0);



        if (d00 >= (1-1e-6)*(1-1e-6))
        {
            if (verbose >= 1)
                System.out.println("DOING THE INFINITESIMAL THING! (no way this will work)");
{ if (!(false)) set_exception("HyperbolicUtils.java","_calcPointOnWeightedSector2" , "VecMath.dot(p0,p0)<1" ); };

            double v[][] = new double[nPoints][];
            int i;
            for (i = 0; (i) < (nPoints); ++i)
                v[i] = _calcEndOfGeodesic2(p0, p[i]);

            double dots[] = new double[nPoints];
            for (i = 0; (i) < (nPoints); ++i)
                dots[i] = VecMath.dot(p0,v[i]);




            double u[][] = new double[nPoints][dim];
            for (i = 0; (i) < (nPoints); ++i)
            {
                VecMath.vpsxv(u[i], v[i], -dots[i]/d00, p0);
                double len = Math.sqrt(VecMath.dot(u[i],u[i]));



                if (len != 0.) VecMath.vxs(u[i],u[i],1./len);
            }

            double w[][] = new double[nPoints][dim];
            for (i = 0; (i) < (nPoints); ++i)
            {
                double H = .5 * Math.sqrt( 2. / (1.-dots[i]) - 1. );
                VecMath.vxs(w[i], u[i], H);
            }



            double c[] = new double[nPoints];
            for (i = 0; (i) < (nPoints); ++i)
            {
                double otherWs[][] = new double[nPoints-1][dim];
                int j;
                for (j = 0; (j) < (nPoints-1); ++j)
                    VecMath.copyvec(otherWs[j],w[(i+j+1)%nPoints]);
                double oppositeSideContent = calcSimplexContent(2, nPoints-1,otherWs);
                c[i] = b[i] * oppositeSideContent;
            }

            double sum = 0.;
            for (i = 0; (i) < (nPoints); ++i)
                sum += c[i];
            { if (!(sum != 0.)) set_exception("HyperbolicUtils.java","_calcPointOnWeightedSector2" , "sum != 0." ); };
            for (i = 0; (i) < (nPoints); ++i)
                c[i] *= 1./sum;


            double w0[] = new double[dim];
            VecMath.vxm(w0, c, w);
            double H = Math.sqrt(VecMath.dot(w0,w0));
            double u0[] = new double[dim];
            if (H == 0.)
                VecMath.zerovec(u0);
            else
                VecMath.vxs(u0, w0, 1./H);



            double x = 2/(4*H*H+1)-1;
            VecMath.sxvpsxv(result, -x, p0,
                                    Math.sqrt(1.-x*x), u0);
            return;
        }


        Isometry2 F, inverseF;
        F = Isometry2.pureTranslation(p0[0],p0[1]);
        inverseF = Isometry2.pureTranslation(-p0[0],-p0[1]);
        double p_[][] = new double[nPoints][dim];
        int i;
        for (i = 0; (i) < (nPoints); ++i)
            inverseF.apply(p[i], p_[i]);


        double c[] = new double[nPoints];
        for (i = 0; (i) < (nPoints); ++i)
        {
            double otherPoints[][] = new double[nPoints-1][];
            int j;
            for (j = 0; (j) < (nPoints-1); ++j)
                otherPoints[j] = p_[(i+j+1)%nPoints];
            double oppositeSideContent = calcContentWithOrigin(dim,nPoints-1,otherPoints);
            c[i] = b[i] * oppositeSideContent;
        }
        VecMath.vxm(result, c, p_);

        double d = VecMath.dot(result, result);


        if (d == 0.)
            d = 1e-20;

        { if (!(d != 0.)) set_exception("HyperbolicUtils.java","_calcPointOnWeightedSector2" , "d != 0." ); };

        VecMath.vxs(result, result, 1./Math.sqrt(d));


        F.apply(result, result);

    }






    public static double
    calcContentWithOrigin(int n, int k, double _v[ ][ ])
    {

        double v[][] = new double[k][n];
        VecMath.copymat(v, _v);



        double prod = 1.;
        int i;
        for (i = 0; (i) < (k); ++i)
        {
            double dotii = VecMath.dot(v[i],v[i]);
            if (dotii == 0.)
                return 0.;
            int j;
            for (j = i+1; j < k; ++j)
            {


                VecMath.vpsxv(v[j], v[j],
                                    -VecMath.dot(v[j],v[i]) * (1./dotii),
                                    v[i]);
            }
            prod *= dotii;
        }
        return Math.sqrt(prod);
    }



    public static double
    calcSimplexContent(int n, int k, double p[ ][ ])
    {
        double v[][] = new double[k-1][n];

        int i;
        for (i = 0; (i) < (k-1); ++i)
            VecMath.vmv(v[i], p[i+1], p[0]);
        return calcContentWithOrigin(n, k-1, v);
    }

    public static boolean pointIsToRightSideOfLinePoincare(Complex p,
                                                           Complex l0,
                                                           Complex l1,
                                                           double eps)
    {




        Isometry2 I = Isometry2.pureTranslation(-l0.x, -l0.y);
        p = I.apply(p);

        l1 = I.apply(l1);

        double twiceTriangleArea = p.x*l1.y - p.y*l1.x;
        return twiceTriangleArea > eps*eps;
    }

    public static boolean pointIsToRightSideOfLinePoincare(double p[ ],
                                                           double l0[ ],
                                                           double l1[ ],
                                                           double eps)
    {
        return pointIsToRightSideOfLinePoincare(new Complex(p),
                                                new Complex(l0),
                                                new Complex(l1),
                                                eps);
    }



    public static void
    calcFaceCenter(double p, int q, double result[ ])
    {


        double s = Math.sin(Math.PI/q);
        double c = Math.cos(Math.PI/q);
        double r = (p==0 ? 1. : h2eNorm(MyMath.acosh(c/s/Math.tan(Math.PI/p))));
        result[0] = r * c;
        result[1] = r * s;
    }

    public static void
    calcFaceCenterFromHalfEdgeLength(double p, double halfEdgeLength, double result[ ])
    {
        double polygonHalfAngle = polygonHalfAngle(p, halfEdgeLength);
        double circumRadius = polygonCircumRadius(p, halfEdgeLength);
        double r = h2eNorm(circumRadius);
        result[0] = r * Math.cos(polygonHalfAngle);
        result[1] = r * Math.sin(polygonHalfAngle);
    }

    public static double
    polygonHalfAngle(double p, double halfEdgeLength)
    {
        return Math.asin((p == 0. ? 1. : Math.cos(Math.PI/p))
                       / MyMath.cosh(halfEdgeLength));
    }

    public static double
    polygonCircumRadius(double p, double halfEdgeLength)
    {
        if (p == 0.)
            return Double.POSITIVE_INFINITY;
        else
            return MyMath.asinh( MyMath.sinh(halfEdgeLength) / Math.sin(Math.PI/p));
    }

    public static double
    polygonInRadius(double p, double halfEdgeLength)
    {
        if (p == 0.)
            return Double.POSITIVE_INFINITY;
        else
            return MyMath.asinh( MyMath.tanh(halfEdgeLength) / Math.tan(Math.PI/p));
    }

    public static double
    calcUniformTilingHalfEdgeLength(double p[], double q, double density)
    {
        int nps = p.length;
        if (nps == 0 || q == 0.)
            return 0.;





        double max_cos_pi_over_p = Double.NEGATIVE_INFINITY;
        for (int i = 0; i < nps; ++i)
        {
            if (p[i] == 0.)
            {
                max_cos_pi_over_p = 1.;
                break;
            }
            double cos_pi_over_p = Math.cos(Math.PI/p[i]);
            if (cos_pi_over_p > max_cos_pi_over_p)
                max_cos_pi_over_p = cos_pi_over_p;
        }

        double lo = 1.;
        double hi = max_cos_pi_over_p/Math.sin(density*Math.PI/(nps*q));
        if (Double.isInfinite(hi))
            hi = 1e9;
        double mid;
        while ((mid = lo + (hi-lo)*.5) != lo && mid != hi)
        {
            double sum = 0.;
            for (int i = 0; i < nps; ++i)
            {

                if (1./(p[i] == 0. ? p[i] : Math.sin(Math.PI/p[i])) < 0.)
                    sum -= q*2.*Math.asin((p[i]==0.?1.:Math.cos(Math.PI/p[i]))/mid);
                else
                    sum += q*2.*Math.asin((p[i]==0.?1.:Math.cos(Math.PI/p[i]))/mid);
            }
            if (sum > density*2*Math.PI)
                lo = mid;
            else
                hi = mid;
        }
        double result = MyMath.acosh(mid);
        { if (!(!Double.isNaN(result))) set_exception("HyperbolicUtils.java","calcUniformTilingHalfEdgeLength" , "!Double.isNaN(result)" ); };
        { if (!(!Double.isInfinite(result))) set_exception("HyperbolicUtils.java","calcUniformTilingHalfEdgeLength" , "!Double.isInfinite(result)" ); };
        return result;
    }
    public static double
    calcUniformTilingHalfEdgeLength(int p[], double q, double density)
    {
        double pdouble[] = new double[p.length];
        for (int i = 0; (i) < (p.length); ++i)
            pdouble[i] = p[i];
        return calcUniformTilingHalfEdgeLength(pdouble, q, density);
    }







    public static double
    calcRegularTilingHalfEdgeLength(double p, int q)
    {

        double cos_pi_p = (p==0 ? 1. : Math.cos(Math.PI/p));
        double sin_pi_q = (q==0 ? 0. : Math.sin(Math.PI/q));
        double result = MyMath.acosh( cos_pi_p / sin_pi_q);
        { if (!(!Double.isNaN(result))) set_exception("HyperbolicUtils.java","calcRegularTilingHalfEdgeLength" , "!Double.isNaN(result)" ); };
        { if (!(!Double.isInfinite(result))) set_exception("HyperbolicUtils.java","calcRegularTilingHalfEdgeLength" , "!Double.isInfinite(result)" ); };
        { if (!(result > 1e-6)) set_exception("HyperbolicUtils.java","calcRegularTilingHalfEdgeLength" , "result > 1e-6" ); };
        return result;
    }






    public static double
    solveTriangleAFromAlphaCBeta(double alpha, double C, double beta)
    {
        double sinAlpha = Math.sin(alpha);
        double temp = MyMath.cosh(C)*sinAlpha*Math.sin(beta) - Math.cos(alpha)*Math.cos(beta);
        double A = MyMath.asinh(MyMath.sinh(C)*sinAlpha / Math.sqrt(1. - temp*temp));
        return A;
    }






    public static double
    solveTriangleCFromAGammaB(double A, double gamma, double B)
    {
        return MyMath.acosh(MyMath.cosh(A)*MyMath.cosh(B)
                          - Math.cos(gamma)*MyMath.sinh(A)*MyMath.sinh(B));
    }






    public static double
    solveTriangleAlphaFromABC(double A, double B, double C)
    {
        return Math.acos((MyMath.cosh(B)*MyMath.cosh(C) - MyMath.cosh(A))
                       / (MyMath.sinh(B)*MyMath.sinh(C)));
    }

    public static double hangle(double p0[ ], double p1[ ], double p2[ ])
    {
        Isometry2 take_p1_to_origin = Isometry2.pureTranslation(-p1[0], -p1[1]);
        double q0[] = take_p1_to_origin.apply(p0);
        double q2[] = take_p1_to_origin.apply(p2);
        VecMath.normalize(q0,q0);
        VecMath.normalize(q2,q2);
        return VecMath.angleBetweenUnitVectors(q0,q2);
    }


    public static boolean findBaryCoordsThatMakeSnubUniform(
             double result[ ],
             double verts[ ][ ],
             int verbose)
    {



        double schwarzNumbers[] = new double[verts.length];
        for (int i = 0; (i) < (verts.length); ++i)
            schwarzNumbers[i] = Math.PI/hangle(verts[(((i-1)+(verts.length))%(verts.length))],
                                               verts[i],
                                               verts[(i+1)%verts.length]);
        double snubFaceList[] = new double[2*verts.length];
        for (int i = 0; (i) < (verts.length); ++i)
        {
            snubFaceList[2*i+0] = verts.length;
            snubFaceList[2*i+1] = schwarzNumbers[i];
        }
        double snubHalfEdgeLength = calcUniformTilingHalfEdgeLength(snubFaceList, 1., 1.);

        double distsToVerts[] = new double[verts.length];
        for (int i = 0; (i) < (verts.length); ++i)
            distsToVerts[i] = polygonCircumRadius(schwarzNumbers[i], snubHalfEdgeLength);

        double resultPoint[] = new double[2];
        if (!hBary2FromDistances(verts.length,
                                 resultPoint,
                                 distsToVerts,
                                 verts))
            return false;
        return unHBary2(verts.length,
                        result,
                        resultPoint,
                        verts,
                        verbose);
    }

    public static boolean findBaryCoordsThatMakeSnubUniform(
             double result[ ],
             Complex complexVerts[ ],
             int verbose)
    {
        double verts[][] = new double[complexVerts.length][];
        for (int i = 0; (i) < (verts.length); ++i)
            verts[i] = new double[]{complexVerts[i].x,complexVerts[i].y};
        return findBaryCoordsThatMakeSnubUniform(
                   result,
                   verts,
                   verbose);
    }



    private static class Node
    {
        public Node(int valence, int sign, int start)
        {
            this.neighbors = Arrays.fill(valence, -1);
            this.backEdgeInds = Arrays.fill(valence, -1);
            this.sign = sign;
            this.start = start;
        }
        int neighbors[ ];
        int backEdgeInds[ ];
        int sign;
        int start;
    }

    public static int
    UnCachedEnumerateIsometry2GroupForUniformTiling(
                            Isometry2 F0,
                            int nGens,
                            Isometry2 gens[ ],
                            int pp[],
                            int rotationalMultiplicity,
                            int backEdgeInds[ ],
                            int maxIsometries,
                            Isometry2 return_isometries[],
                            int return_isometry_composition_lengths[],
                            Isometry2 return_smallestIsometry,
                            int maxLevels,
                            double tooFar,
                            double tooSmall,
                            int verbose)
    {
        long startTime = System.currentTimeMillis();

        { if (!(rotationalMultiplicity >= 1)) set_exception("HyperbolicUtils.java","UnCachedEnumerateIsometry2GroupForUniformTiling" , "rotationalMultiplicity >= 1" ); };

        double smallestTranslateSqrd = (((F0.P.x)*(F0.P.x))+((F0.P.y)*(F0.P.y)));
        return_smallestIsometry.set(F0);

        if (maxLevels == 0) return 0;
        if (maxIsometries == 0) return 0;



        int backEdgeSigns[] = new int[backEdgeInds.length];
        {
            backEdgeInds = (int[])Arrays.copy(backEdgeInds, 1);
            for (int i = 0; (i) < (backEdgeSigns.length); ++i)
            {
                backEdgeSigns[i] = (backEdgeInds[i] < 0 ? -1 : 1);
                if (backEdgeSigns[i] == -1)
                    backEdgeInds[i] ^= ~0;

                backEdgeSigns[i] *= -1;
            }
        }



        {
            int expandedPP[] = new int[pp.length * rotationalMultiplicity];
            int expandedBackEdgeInds[] = new int[pp.length * rotationalMultiplicity];
            int expandedBackEdgeSigns[] = new int[pp.length * rotationalMultiplicity];
            for (int i = 0; (i) < (expandedPP.length); ++i)
            {
                expandedPP[i] = pp[i%pp.length];
                expandedBackEdgeSigns[i] = backEdgeSigns[i%pp.length];
                expandedBackEdgeInds[i] = backEdgeInds[i%pp.length] + i/pp.length*pp.length;
            }

            pp = expandedPP;
            backEdgeInds = expandedBackEdgeInds;
            backEdgeSigns = expandedBackEdgeSigns;
            rotationalMultiplicity = 1;
        }

        { if (!(pp.length == nGens)) set_exception("HyperbolicUtils.java","UnCachedEnumerateIsometry2GroupForUniformTiling" , "pp.length == nGens" ); };

        Node nodes[] = new Node[maxIsometries];
        int n = 0;
        if (n >= maxIsometries) return n;
        return_isometries[n] = F0;
        nodes[n] = new Node(nGens, 1, 0);
        if (return_isometry_composition_lengths != null)
            return_isometry_composition_lengths[0] = 0;
        n++;
        if (n >= maxIsometries) return n;

        int levelStart = 0;
        int iNode = 0;
        int iLevel = 0;
        for (iLevel = 0; (iLevel) < (maxLevels); ++iLevel)
        {



            int levelEnd = n;
            for (; iNode < levelEnd; ++iNode)
            {




                for (int _iGen = 0; (_iGen) < (nGens); ++_iGen)
                {
                    int iGen = (((nodes[iNode].start + nodes[iNode].sign*_iGen)+(nGens))%(nGens));

                    if (verbose >= 3)
                        System.out.println("    Trying "+iNode+":"+iGen+"");
                    if (nodes[iNode].neighbors[iGen] != -1)
                    {
                        if (verbose >= 3)
                            System.out.println("        (ha, "+iNode+":"+iGen+" -> "+nodes[iNode].neighbors[iGen]+":"+nodes[iNode].backEdgeInds[iGen]+")");
                        continue;
                    }





                    {
                        for (int dir = 1; dir >= -1; dir -= 2)
                        {
                            int jNode = iNode;
                            int jGen = (((iGen + dir*nodes[jNode].sign)+(nGens))%(nGens));
                            int faceSize = pp[dir*nodes[iNode].sign==-1 ? (((iGen-1)+(nGens))%(nGens)) : iGen];
                            if (verbose >= 3)
                                System.out.print("            ");
                            for (int iAlongFace = 0; (iAlongFace) < (faceSize-1); ++iAlongFace)
                            {

                                Node node = nodes[jNode];
                                jNode = node.neighbors[jGen];
                                jGen = node.backEdgeInds[jGen];
                                if (verbose >= 3)
                                    System.out.print(""+jNode+":"+jGen+" ");
                                if (jNode == -1)
                                    break;
                                jGen = (((jGen + dir*nodes[jNode].sign)+(nGens))%(nGens));
                                if (verbose >= 3)
                                    System.out.print(""+jNode+":"+jGen+" ");
                            }
                            if (verbose >= 3)
                                System.out.println();
                            if (jNode != -1)
                            {
                                if (verbose >= 3)
                                    System.out.println("        (HA! found it!!! "+iNode+":"+iGen+" -> "+jNode+":"+jGen+" (dir="+dir+"))");
                                nodes[iNode].neighbors[iGen] = jNode;
                                nodes[iNode].backEdgeInds[iGen] = jGen;
                                nodes[jNode].neighbors[jGen] = iNode;
                                nodes[jNode].backEdgeInds[jGen] = iGen;
                                break;
                            }
                        }
                        if (nodes[iNode].neighbors[iGen] != -1)
                            continue;
                    }


                    Isometry2 I = return_isometries[iNode];
                    Isometry2 F = Isometry2.mul(I, gens[iGen], return_isometries[n]);


                    if ((gens[iGen].P.x != 0. || gens[iGen].P.y != 0.)
                     && (((F.P.x-I.P.x)*(F.P.x-I.P.x))+((F.P.y-I.P.y)*(F.P.y-I.P.y)))
                                               < tooSmall*tooSmall)
                    {
                        if (verbose >= 2)
                        {
                            System.out.print("S");
                        }
                        continue;
                    }


                    double thisTranslateSqrd = (((F.P.x)*(F.P.x))+((F.P.y)*(F.P.y)));
                    if (thisTranslateSqrd >= tooFar*tooFar)
                    {
                        if (verbose >= 2)
                        {
                            System.out.print("F");
                        }
                        continue;
                    }

                    int jNode = n;
                    int jGen = backEdgeInds[iGen];
                    nodes[jNode] = new Node(nGens,
                                        nodes[iNode].sign * backEdgeSigns[iGen],
                                        jGen);
                    nodes[iNode].neighbors[iGen] = jNode;
                    nodes[iNode].backEdgeInds[iGen] = jGen;
                    nodes[jNode].neighbors[jGen] = iNode;
                    nodes[jNode].backEdgeInds[jGen] = iGen;
                    if (return_isometry_composition_lengths != null)
                        return_isometry_composition_lengths[jNode] = return_isometry_composition_lengths[iNode]+1;

                    if (return_isometry_composition_lengths != null
                        ? return_isometry_composition_lengths[jNode]%2 == 0
                        : return_isometries[jNode].R == 1)
                    {
                        if (thisTranslateSqrd < smallestTranslateSqrd)
                        {
                            smallestTranslateSqrd = thisTranslateSqrd;
                            return_smallestIsometry.set(F);
                        }
                    }

                    if (verbose >= 3)
                        System.out.println("        adding "+iNode+" -> "+nodes[iNode].neighbors[iGen]+"");

                    n++;

                    if (n >= maxIsometries) break;
                }
                if (n >= maxIsometries) break;
            }
            if (n >= maxIsometries) break;

            levelStart = levelEnd;
        }

        long endTime = System.currentTimeMillis();

        if (verbose >= 1)
        {

            if (iNode == n)
                System.out.println("Got all "+n+" return_isometries\n");
            if (iLevel >= maxLevels)
                System.out.println("Terminating enumeration at "+iLevel+" levels\n");
            if (n == maxIsometries)
                System.out.println("Terminating enumeration at "+n+" return_isometries\n");
            {System.out.print("    ");System.out.println("endTime-startTime" + " = " + (endTime-startTime));};
        }

        return n;
    }




    public static int
    UnCachedEnumerateIsometry2Group(
                            Isometry2 F0,
                            int nGens,
                            Isometry2 gens[ ],
                            int rotationalMultiplicity,
                            int maxIsometries,
                            Isometry2 return_isometries[],
                            Isometry2 return_smallestIsometry,
                            int maxLevels,
                            double tooFar,
                            double tooSmall,
                            int verbose)
    {
        long startTime = System.currentTimeMillis();

        { if (!(rotationalMultiplicity >= 1)) set_exception("HyperbolicUtils.java","UnCachedEnumerateIsometry2Group" , "rotationalMultiplicity >= 1" ); };
        Isometry2HashTable hashTable = new Isometry2HashTable();


        double smallestTranslateSqrd = (((F0.P.x)*(F0.P.x))+((F0.P.y)*(F0.P.y)));
        return_smallestIsometry.set(F0);

        if (maxLevels == 0) return 0;
        if (maxIsometries == 0) return 0;

        int n = 0;
        if (n >= maxIsometries) return n;
        boolean inserted = hashTable.Insert(F0);
        { if (!(inserted)) set_exception("HyperbolicUtils.java","UnCachedEnumerateIsometry2Group" , "inserted" ); };
        return_isometries[n++] = F0;
        if (n >= maxIsometries) return n;
        int i;

        int nIsometriesOnPreviousLevels = 1;
        int generatingLevel = 1;
        if (generatingLevel >= maxLevels)
            return n;

        Isometry2 rotation = Isometry2.pureRotation(Math.PI/rotationalMultiplicity);

        Isometry2 temp = new Isometry2();

if (verbose >= 1) {System.out.print("    ");System.out.println("rotationalMultiplicity" + " = " + (rotationalMultiplicity));};
rotationalMultiplicity = 1;
if (verbose >= 1) {System.out.print("    ");System.out.println("rotationalMultiplicity" + " = " + (rotationalMultiplicity));};

        for (i = 0; (i) < (n); ++i)
        {
            if (i == nIsometriesOnPreviousLevels)
            {
                if (++generatingLevel >= maxLevels)
                {
                    break;
                }
                nIsometriesOnPreviousLevels = n;
            }
            Isometry2 I = new Isometry2(return_isometries[i]);
            int iRot;
            for (iRot = 0; (iRot) < (rotationalMultiplicity); ++iRot)
            {
                int iGen;
                for (iGen = 0; (iGen) < (nGens); ++iGen)
                {

                    Isometry2 F = Isometry2.mul(I, gens[iGen], return_isometries[n]);


                    if ((gens[iGen].P.x != 0. || gens[iGen].P.y != 0.)
                     && (((F.P.x-I.P.x)*(F.P.x-I.P.x))+((F.P.y-I.P.y)*(F.P.y-I.P.y)))
                                               < tooSmall*tooSmall)
                    {
                        if (verbose >= 2)
                        {
                            System.out.print("S");
                        }
                        continue;
                    }


                    double thisTranslateSqrd = (((F.P.x)*(F.P.x))+((F.P.y)*(F.P.y)));
                    if (thisTranslateSqrd >= tooFar*tooFar)
                    {
                        if (verbose >= 2)
                        {
                            System.out.print("F");
                        }
                        continue;
                    }

                    inserted = hashTable.Insert(F);
                    if (!inserted)
                        continue;




                    if (thisTranslateSqrd < smallestTranslateSqrd)
                    {
                        smallestTranslateSqrd = thisTranslateSqrd;
                        return_smallestIsometry.set(F);
                    }

                    n++;

                    if (n >= maxIsometries) break;
                }

                if (n >= maxIsometries) break;

                if (iRot == rotationalMultiplicity-1)
                    break;


                Isometry2.mul(I, rotation, temp);
                I.set(temp);
            }
            if (n >= maxIsometries) break;
        }

        if (verbose >= 1)
        {

            if (i == n)
                System.out.println("Got all "+n+" return_isometries\n");
            if (generatingLevel >= maxLevels)
                System.out.println("Terminating enumeration at "+generatingLevel+" levels\n");
            if (n == maxIsometries)
                System.out.println("Terminating enumeration at "+n+" return_isometries\n");
        }

        long endTime = System.currentTimeMillis();
        if (verbose >= 1)
        {
            {System.out.print("    ");System.out.println("endTime-startTime" + " = " + (endTime-startTime));};
            {System.out.print("    ");System.out.println("hashTable.nInsertSuccesses" + " = " + (hashTable.nInsertSuccesses));};
            {System.out.print("    ");System.out.println("hashTable.nInsertFailures" + " = " + (hashTable.nInsertFailures));};
        }

        return n;
    }

    public static int
    subdivideArc(double p0[ ], double p1[ ], double tolerance, int maxIntermediatePoints, int startI, double result_intermediatePoints[][ ])
    {
        double fullDist;
        {


            Complex _p0 = new Complex(p0);
            Complex _p1 = new Complex(p1);
            Complex _p1_minus_p0 = Isometry2.pureTranslation(-_p0.x,-_p0.y).apply(_p1);
            fullDist = e2hNorm(Math.sqrt((((_p1_minus_p0.x)*(_p1_minus_p0.x))+((_p1_minus_p0.y)*(_p1_minus_p0.y)))));
        }


        final int N = p0.length;
        { if (!(N == 2)) set_exception("HyperbolicUtils.java","UnCachedEnumerateIsometry2Group" , "N == 2" ); };


        int nIntermediatePoints = maxIntermediatePoints;
        int i;
        double sinh_fullDist = MyMath.sinh(fullDist);
        for (i = 0; (i) < (nIntermediatePoints); ++i)
        {
            double t = (double)(i+1)/(nIntermediatePoints+1);
            if (false)
                VecMath.sxvpsxv(result_intermediatePoints[startI+i],
                                MyMath.sinh((1-t)*fullDist)/sinh_fullDist, p0,
                                MyMath.sinh( t*fullDist)/sinh_fullDist, p1);
            else
                hlerp2_poincare(p0, p1, t, result_intermediatePoints[startI+i]);
        }


        return nIntermediatePoints;
    }



    static double hlerp2_poincare(double p0[], double p1[], double t, double result[])
    {
        Isometry2 take_p0_to_origin = Isometry2.pureTranslation(-p0[0],-p0[1]);
        Isometry2 take_origin_to_p0 = Isometry2.pureTranslation(p0[0],p0[1]);
        double p1_[] = take_p0_to_origin.apply(p1);
        double normSqrd = (((p1_[0])*(p1_[0]))+((p1_[1])*(p1_[1])));
        if (normSqrd < 1e-12)
        {
            if (result != null)
            {
                result[0] = p0[0] + t*(p1[0]-p0[0]);
                result[1] = p0[1] + t*(p1[1]-p0[1]);
            }
            return 2*Math.sqrt(normSqrd);
        }
        if (normSqrd > 1.)
        {
            if (result != null)
            {
                result[0] = 0.;
                result[1] = 0.;
            }
            return -1.;
        }
        double norm = Math.sqrt(normSqrd);
        double atanNorm = MyMath.atanh(norm);
        if (result != null)
        {
            double desiredNorm = MyMath.tanh(t * atanNorm);

            double p_[] = p1_;
            VecMath.vxs(p_, p1_, desiredNorm/norm);
            take_origin_to_p0.apply(p_, result);
        }
        return 2*atanNorm;
    }


    static double hlerp2_poincare(Complex p0, Complex p1, double t, Complex result)
    {
        double temp[] = new double[2];
        double dist = hlerp2_poincare(new double[]{p0.x,p0.y}, new double[]{p1.x,p1.y}, t, temp);
        result.set(temp[0],temp[1]);
        return dist;
    }

    static double hdist2_poincare(double p0[ ], double p1[ ])
    {
        return hlerp2_poincare(p0, p1, 0., null);
    }



    static double hlerp2_poincare(double p0[], double p1[], double t, double aspect, double result[])
    {
        if (aspect == 1.)
            return hlerp2_poincare(p0,p1,t,result);

        double _p0[] = new double[2];
        double _p1[] = new double[2];
        VecMath.vxs(_p0, p0, Math.pow(Math.sqrt(VecMath.dot(p0,p0)), 1./aspect-1.));
        VecMath.vxs(_p1, p1, Math.pow(Math.sqrt(VecMath.dot(p1,p1)), 1./aspect-1.));
        double dist = hlerp2_poincare(_p0,_p1,t,result);
        if (result != null)
            VecMath.vxs(result, result, Math.pow(Math.sqrt(VecMath.dot(result,result)), aspect-1.));
        return dist;
    }


    static double hlerp2_uhp(double p0[], double p1[], double t,
                              double aspect,
                              double result[])
    {
        double p0_[] = new double[2];
        double p1_[] = new double[2];
        double p_[] = (result == null ? null : new double[2]);


        uhp2d(p0[0], p0[1]/aspect, p0_);
        uhp2d(p1[0], p1[1]/aspect, p1_);


        double ret = hlerp2_poincare(p0_, p1_, t, p_);

        if (result != null)
        {

            d2uhp(p_[0], p_[1], result);

            result[1] *= aspect;
        }

        return ret;
    }

    static double hlerp2_simple_inversion(double p0[], double p1[], double t,
                                double aspect,
                                boolean conformal,
                                double result[])
    {


        double p0_[] = new double[2];
        double p1_[] = new double[2];
        double u0[] = new double[2];
        {
            double invLength0 = 1./Math.sqrt(VecMath.dot(p0,p0));
            double invLength1 = 1./Math.sqrt(VecMath.dot(p1,p1));
            double u1[] = new double[2];
            VecMath.vxs(u0, p0, invLength0);
            VecMath.vxs(u1, p1, invLength1);
            double u[] = new double[2];
            VecMath.vpv(u, u0, u1);
            VecMath.vxs(u, u, 1./Math.sqrt(VecMath.dot(u,u)));
            VecMath.vxs(p1_, u, Math.sqrt(invLength1));
            VecMath.vxs(p0_, u0, Math.sqrt(invLength0));
        }



        double ret = hlerp2_poincare(p0_, p1_, t, result);

        if (result != null && ret != -1.)
        {
            double finalLength = 1./VecMath.dot(result,result);
            double u[] = new double[2];
            VecMath.vxs(u, result, Math.sqrt(finalLength));
            double halfway[] = new double[2];
            VecMath.vxs(halfway, u, VecMath.dot(u0,u));
            double fullway[] = new double[2];
            VecMath.sxvpsxv(fullway, 2., halfway, -1., u0);
            VecMath.vxs(result, fullway, finalLength);
        }
        return ret;
    }



    static double
    distToCircleCenteredAtMinus1(double x, double y)
    {

        return (x*x + y*(y+2.)) / (Math.sqrt(x*x + (y+1.)*(y+1.)) + 1.);
    }


    static double hlerp2_polar(double p0[], double p1[], double t,
                               double aspect,
                               boolean conformal,
                               double result[])
    {



        double x0 = p0[0];
        double y0 = p0[1]-1.;
        double x1 = p1[0];
        double y1 = p1[1]-1.;
        double actualAlt0 = distToCircleCenteredAtMinus1(x0,y0);
        double actualAlt1 = distToCircleCenteredAtMinus1(x1,y1);







        double angle0 = Math.atan2(-x0,y0+1.);
        double angle1 = Math.atan2(-x1,y1+1.);



        {
            if (angle1 > angle0 + Math.PI)
                angle1 -= 2*Math.PI;
            else if (angle1 < angle0 - Math.PI)
                angle1 += 2*Math.PI;
        }

        double alt0 = MyMath.log1p(actualAlt0);
        double alt1 = MyMath.log1p(actualAlt1);
        if (!conformal)
        {
            alt0 = actualAlt0;
            alt1 = actualAlt1;
        }

        double p0_[] = {angle0, alt0};
        double p1_[] = {angle1, alt1};
        double ret = hlerp2_uhp(p0_, p1_, t, aspect, result);
        {
            if (result != null)
            {
                double angle = result[0];
                double alt = result[1];
                double actualAlt = conformal ? MyMath.expm1(alt)
                                             : alt;
                double r = actualAlt+1.;




                double x = -r * Math.sin(angle);


                double y = actualAlt - r * MyMath.cosf1(angle);

                result[0] = x;
                result[1] = y+1.;
            }
        }

        return ret;
    }

    static boolean hspline(double t0, double p0[],
                           double t1, double p1[],
                           double t2, double p2[],
                           double t3, double p3[],
                           double t,
                           double aspect,
                           boolean conformal,
                           LerpFunction lerp,
                           double result[])
    {
        final int verbose = 0;
        if (verbose >= 1)
        {
            System.out.println();
            System.out.println("in hspline");
        }




        {
            double transT = -t1;
            double scaleT = 1./(t2-t1);
            t0 = (t0+transT)*scaleT;
            t1 = 0.;
            t2 = 1.;
            t3 = (t3+transT)*scaleT;
            t = (t+transT)*scaleT;
        }




        double W[][] = new double[4][4];
        {
            double invT[][] = {
                {1,0,-3,2},
                {0,0,3,-2},
                {0,1,-2,1},
                {0,0,-1,1},
            };
            double thatHairyMatrixOnTheRight[][] = {
                {0, 0, -1/((1-t0)*(-t0)), 0 },
                {1, 0, 1/((1-t0)*(-t0))-(-t0)/(1-t0), -(t3-1)/t3 },
                {0, 1, (-t0)/(1-t0), (t3-1)/t3-1/(t3*(t3-1))},
                {0, 0, 0, 1/(t3*(t3-1)) },
            };
            if (verbose >= 2) System.out.println("    " + "thatHairyMatrixOnTheRight" + " =\n" + VecMath.toString(thatHairyMatrixOnTheRight));
            VecMath.mxm(W, thatHairyMatrixOnTheRight, invT);
            if (verbose >= 2) System.out.println("    " + "W" + " =\n" + VecMath.toString(W));

            if (verbose >= 2)
            {
                double T[][] = {
                    {1,1,0,0},
                    {0,1,1,1},
                    {0,1,0,2},
                    {0,1,0,3},
                };
                double WT[][] = new double[4][4];
                VecMath.mxm(WT, W, T);
                System.out.println("should be H...");
                System.out.println("    " + "WT" + " =\n" + VecMath.toString(WT));

                double T_invT[][] = new double[4][4];
                double invT_T[][] = new double[4][4];
                VecMath.mxm(T_invT, T, invT);
                VecMath.mxm(invT_T, invT, T);
                System.out.println("    " + "T_invT" + " =\n" + VecMath.toString(T_invT));
                System.out.println("    " + "invT_T" + " =\n" + VecMath.toString(invT_T));
            }
        }





        double w0, w1, w2, w3;
        {



            double tVec[] = {1, t, t*t, t*t*t};
            double wVec[] = new double[4];
            VecMath.mxv(wVec, W, tVec);
            w0 = wVec[0];
            w1 = wVec[1];
            w2 = wVec[2];
            w3 = wVec[3];

        }

        {
            double w1_01 = (w0+w1 == 0. ? .5 : w1/(w0+w1));
            double w3_23 = (w2+w3 == 0. ? .5 : w3/(w2+w3));
            double w23 = w2+w3;
            double A[] = new double[2];
            double B[] = new double[2];

            if (verbose >= 1)
            {
                System.out.println("    " + "p0" + " = " + VecMath.toString(p0));
                System.out.println("    " + "p1" + " = " + VecMath.toString(p1));
                System.out.println("    " + "p2" + " = " + VecMath.toString(p2));
                System.out.println("    " + "p3" + " = " + VecMath.toString(p3));
                {System.out.print("    ");System.out.println("t" + " = " + (t));};
                {System.out.print("    ");System.out.println("w0" + " = " + (w0));};
                {System.out.print("    ");System.out.println("w1" + " = " + (w1));};
                {System.out.print("    ");System.out.println("w2" + " = " + (w2));};
                {System.out.print("    ");System.out.println("w3" + " = " + (w3));};
                {System.out.print("    ");System.out.println("w1_01" + " = " + (w1_01));};
                {System.out.print("    ");System.out.println("w3_23" + " = " + (w3_23));};
                {System.out.print("    ");System.out.println("w23" + " = " + (w23));};
            }

            if (lerp.apply(p0,p1, w1_01, aspect, conformal, A) == -1.) return false;
            if (verbose >= 1) System.out.println("    " + "A" + " = " + VecMath.toString(A));
            if (lerp.apply(p2,p3, w3_23, aspect, conformal, B) == -1.) return false;
            if (verbose >= 1) System.out.println("    " + "B" + " = " + VecMath.toString(B));
            if (lerp.apply(A,B, w23, aspect, conformal, result) == -1.) return false;
            if (verbose >= 1) System.out.println("    " + "result" + " = " + VecMath.toString(result));
        }
        return true;
    }



    private static void canonicalCircularInvert(double x, double y, double result[ ])
    {
        double scale = 2. / (x*x + (y+1.)*(y+1.));
        result[0] = x * scale;
        result[1] = (y+1.) * scale - 1.;
    }

    static boolean hlerp3_polar(double p0[], double p1[], double t,
                                double aspect,
                                double result)
    {

        return false;

    }


    static void d2uhp(double x, double y, double result[ ])
    {
        canonicalCircularInvert(x, -y, result);
    }
    static void uhp2d(double x, double y, double result[ ])
    {
        canonicalCircularInvert(x, y, result);
        result[1] *= -1;
    }


    public static final int POINCARE_DISK = 0;
    public static final int KLEIN_DISK = 1;
    public static final int UPPER_HALF_PLANE = 2;
    public static final int ANTIKLEIN_DISK = 3;

     public abstract static class LerpFunction
     {




         public abstract double apply(
                double p0[], double p1[], double t,
                double aspect,
                boolean conformal,
                double result[]);





        public double dist(double p0[], double p1[], double aspect, boolean conformal)
        {
            return apply(p0, p1, .5, aspect, conformal, null);
        }
     }
     public static class Hlerp2_poincare extends LerpFunction
     {
        public double apply(double p0[], double p1[], double t,
                            double aspect,
                            boolean conformal,
                            double result[])
        { return hlerp2_poincare(p0,p1,t,aspect,result); }
     }
     public static class Hlerp2_uhp extends LerpFunction
     {
        public double apply(double p0[], double p1[], double t,
                            double aspect,
                            boolean conformal,
                            double result[])
        { return hlerp2_uhp(p0,p1,t,aspect,result); }
     }
     public static class Hlerp2_polar extends LerpFunction
     {
        public double apply(double p0[], double p1[], double t,
                            double aspect,
                            boolean conformal,
                            double result[])
        { return hlerp2_polar(p0,p1,t,aspect,conformal,result); }
     }
     public static class Hlerp2_simple_inversion extends LerpFunction
     {
        public double apply(double p0[], double p1[], double t,
                            double aspect,
                            boolean conformal,
                            double result[])
        { return hlerp2_simple_inversion(p0,p1,t,aspect,conformal,result); }
     }
     public static class Hlerp2_euclidean extends LerpFunction
     {
        public double apply(double p0[], double p1[], double t,
                            double aspect,
                            boolean conformal,
                            double result[])
        {
            if (result != null)
                VecMath.sxvpsxv(result, 1.-t, p0, t, p1);
            return Math.sqrt(VecMath.distsqrd(p0,p1));
        }
     }
}
