// WeierstrassP
// 
// A Java class to iterate Weierstrass elliptic P functions.
// By Mark McClure
// April, 2008
//
// Primarily meant to be called by Mathematica as a JLink program.
// 
// The algorithm behind wp that computes values of the Weierstrass
// elliptic functions comes from:
// "Iterative Method for Calculation of the Weierstrass Elliptic Function" 
// by R. Coquereaux, A. Grossmann and B. E. Lautrup 
// IMA Journal of Numerical Analysis 1990 10(1):119-128.
//
// The other functions are pretty straightforward techniques in 
// complex dynamics.

// Just learned a couple of things on numerical Java from 
// http://www2.informatik.uni-erlangen.de/Forschung/Publikationen/download/cise-ron.pdf
//  * JVM implementation very important
//  * do for(int i; ...)  Creates a local variable
//  Compress multi-dimensional arrays (A[i][j][k] -> Aij[k])


// Computations are done in terms of complex numbers defined in 
// the Netlib complex class.
import ORG.netlib.math.complex.*;


// The main class.
public class WeierstrassP { 
    static final int DEFAULT_BAIL = 10;

    // Fake main method.
    public static void main (String[] args) {

        int i,j;
        Complex zMin = new Complex(-2.2,-2.2);
        Complex zMax = new Complex(2.2,2.2);
        Complex g2 = new Complex(5.746, 0.0);
        Complex g3 = new Complex(0.0,0.0);
        Complex[] fps = new Complex[1];
        fps[0] = g2.sqrt().scale(0.5);
        int[][][] result;
        
        result = itCounts(zMin,zMax, g2,g3,10,fps, 100, 0.01, 10, 10);
        
        for(i=0; i<11; i++) {
            for(j=0; j<11; j++) {
                System.out.print(result[i][j][0]);
                System.out.print(" - ");
                System.out.print(result[i][j][1]);
                System.out.print("  ---  ");
            }
            System.out.print("\n");
        }
    }

    // wp computes values of the WeierstrassP function.
    public static Complex wp(Complex z, double g2, double g3, int n) {
        int i;
        Complex z0 = new Complex();
        Complex z0sq = new Complex();
        Complex term1 = new Complex();
        Complex term2 = new Complex();
        Complex term3 = new Complex();
        Complex term4 = new Complex();
        Complex out = new Complex();
        
        z0 = z.scale(Math.pow(2,-n));
        z0sq = Complex.pow(z0,2);
        term1 = z0sq.scale(g3/28.0).add(new Complex(g2/20.0, 0));
        out = Complex.pow(z0, -2).add(z0sq.mul(term1));
        for(i=0; i<n; i++) {
            term2 = Complex.pow(Complex.pow(out,2.0).scale(6.0).sub(new Complex(g2/2.0,0.0)),2.0);
            term3 = Complex.pow(out,3.0).scale(4.0).sub(out.scale(g2)).sub(new Complex(g3,0.0));
            term4 = term2.div(term3.scale(4.0));
            out = term4.sub(out.scale(2.0));
        }
        return out;
    }
    
    public static Complex wp(Complex z, Complex g2, Complex g3, int n) {
        int i;
        Complex z0 = new Complex();
        Complex z0sq = new Complex();
        Complex term1 = new Complex();
        Complex term2 = new Complex();
        Complex term3 = new Complex();
        Complex term4 = new Complex();
        Complex out = new Complex();
        
        z0 = z.scale(Math.pow(2,-n));
        z0sq = Complex.pow(z0,2);
        term1 = z0sq.mul(g3.scale(1/28.0)).add(g2.scale(0.05));
        out = Complex.pow(z0, -2).add(z0sq.mul(term1));
        for(i=0; i<n; i++) {
            term2 = Complex.pow(Complex.pow(out,2.0).scale(6.0).sub(g2.scale(0.5)),2.0);
            term3 = Complex.pow(out,3.0).scale(4.0).sub(out.mul(g2)).sub(g3);
            term4 = term2.div(term3.scale(4.0));
            out = term4.sub(out.scale(2.0));
        }
        return out;
    }

    // The critical values of P(z; g2, g3) are the roots of p(z) = 4*x^3 - g2*x -g3.
    // This function computes those roots using a formula based on the general cubic.
    public static Complex[] wpCriticalValues(Complex g2, Complex g3) {
        // Common terms
        Complex DSqrt = (Complex.pow(g3,2.0).scale(27.0).sub(Complex.pow(g2,3.0))).sqrt();
        Complex cubeRootTerm = Complex.pow((g3.scale(9.0)).add(DSqrt.scale(Math.sqrt(3.0))), 1/3.0);
        
        // The first root
        Complex term1 = g2.div(cubeRootTerm.scale(2*Math.pow(3,1/3.0)));
        Complex term2 = cubeRootTerm.scale(1/(2*Math.pow(3,2/3.0)));
        Complex root1 = term1.add(term2);
        
        // Second root
        term1 = g2.mul(new Complex(-1.0,-Math.sqrt(3.0))).div(cubeRootTerm.scale(4*Math.pow(3,1/3.0)));
        term2 = cubeRootTerm.mul(new Complex(-1.0,Math.sqrt(3.0))).scale(1/(4*Math.pow(3,2/3.0)));
        Complex root2 = term1.add(term2);
        
        // Third and final root
        term1 = g2.mul(new Complex(-1.0,Math.sqrt(3.0))).div(cubeRootTerm.scale(4*Math.pow(3,1/3.0)));
        term2 = cubeRootTerm.mul(new Complex(-1.0,-Math.sqrt(3.0))).scale(1/(4*Math.pow(3,2/3.0)));
        Complex root3 = term1.add(term2);
        
        // Return the result
        Complex[] result = new Complex[3];
        result[0] = root1;
        result[1] = root2;
        result[2] = root3;
        return result;
    }


    // Computes the shortest distance between the complex number z0 and
    // the list of complex numbers zs.  Used to test when an iterate is 
    // close enough to an attractive orbit to assume convergence.
    public static double[] dist(Complex z0, Complex[] zs) {
        int i;
        double m;
        int index;
        double thisDist;
        double[] result = new double[2];
        
        m = z0.sub(zs[0]).abs();
        index = 0;
        for(i=1; i<zs.length; i++) {
            thisDist = z0.sub(zs[i]).abs();
            if(m > thisDist) {
                m = thisDist;
                index = i;
            }
        }
        result[0] = m;
        result[1] = new Double(index);
        return result;
    }
        
    // Iterates the Weierstrass elliptic function p(z; g2, g3) at most
    // bail times starting from z0 until dist(z, fps) < tol (in particular,
    // this assumes that the attractive orbits are already known.  Returns
    // an array [number of iterates, fixed point index].
//     public static int[] itCount(Complex z0, double g2, double g3, int n,
//         Complex[] fps, int bail, double tol) {
//     
//         int i;
//         double[] currentDist;
//         Complex z = new Complex(z0);
//         int[] result = new int[2];
//         
//         i = 0;
//         currentDist = dist(z, fps);
//         while(i++<bail && currentDist[0]>tol) {
//             z = wp(z, g2, g3, n);
//             currentDist = dist(z, fps);
//         }
//         result[0] = i;
//         result[1] = (int) currentDist[1];
//         return result;
//     }
    public static int[] itCount(Complex z0, Complex g2, Complex g3, int n,
        Complex[] fps, int bail, double tol) {
    
        int i;
        double[] currentDist;
        Complex z = new Complex(z0);
        int[] result = new int[2];
        
        i = 0;
        currentDist = dist(z, fps);
        while(i++<bail && currentDist[0]>tol) {
            z = wp(z, g2, g3, n);
            currentDist = dist(z, fps);
        }
        result[0] = i;
        result[1] = (int) currentDist[1];
        return result;
    }

//  This version returns complex values, so I can check the iterates as well.
//     public static Complex[] CriticalITCountSquare(Complex g2, int n,
//         int bail, double tol, int maxOrbitLength) {
//     
//         int i;
//         Complex[] recentOrbit = new Complex[maxOrbitLength];
//         double[] currentDist;
//         Complex zeroComplex = new Complex(0.0,0.0);
//         Complex[] result = new Complex[4];
//         
//         i = 0;
//         Complex z = g2.sqrt().scale(0.5);
//         while(i < maxOrbitLength) {
//             z = wp(z, g2, zeroComplex, n);
//             recentOrbit[i] = z;
//             i++;
//         }
//         
//         z = wp(z, g2, zeroComplex, n);
//         currentDist = dist(z, recentOrbit);
//         recentOrbit[i % maxOrbitLength] = z;
//         int  out1 = i % maxOrbitLength;
//         while(i<bail && currentDist[0]>tol) {
//             z = wp(z, g2, zeroComplex, n);
//             currentDist = dist(z, recentOrbit);
//             i = i+1;
//             recentOrbit[i % maxOrbitLength] = z;
//         }
//         result[0] = new Complex(i,0);
//         if((i % maxOrbitLength) > currentDist[1]) {
//             result[1] = new Complex(i % maxOrbitLength - (int) currentDist[1], 0);
//             result[2] = new Complex(1,0);
//         }
//         if((i % maxOrbitLength) <= currentDist[1]) {
//             result[1] = new Complex(maxOrbitLength - ((int) currentDist[1] - i % maxOrbitLength),0);
//             result[2] = new Complex(2,0);
//         }
//         result[3] = z;
//         return result;
//     }

    
    // Loops through all complex numbers in a rectangular grid and
    // performs the iteration of the previous function.
//     public static int[][][] itCounts(Complex zMin, Complex zMax, 
//         double g2, double g3, int n, Complex[] fps, int bail, double tol,
//         int resRe, int resIm) {
//         
//         Complex z = new Complex();
//         int j, k;
//         int[][][] counts = new int[resIm + 1][resRe + 1][2];
//         
//         for(k=0; k <= resIm; k++) {
//             for(j=0; j <= resRe; j++) {
//                 z = new Complex((double) (zMax.re() - zMin.re())*j/resRe + zMin.re(),
//                     (double) (zMax.im() - zMin.im())*k/resIm + zMin.im());
//                 counts[k][j] = itCount(z, g2, g3, n, fps, bail, tol);
//             }
//         }
//         return counts;
//     }
    public static int[][][] itCounts(Complex zMin, Complex zMax, 
        Complex g2, Complex g3, int n, Complex[] fps, int bail, double tol,
        int resRe, int resIm) {
        
        Complex z = new Complex();
        int j, k;
        int[][][] counts = new int[resIm + 1][resRe + 1][2];
        
        for(k=0; k <= resIm; k++) {
            for(j=0; j <= resRe; j++) {
                z = new Complex((double) (zMax.re() - zMin.re())*j/resRe + zMin.re(),
                    (double) (zMax.im() - zMin.im())*k/resIm + zMin.im());
                counts[k][j] = itCount(z, g2, g3, n, fps, bail, tol);
            }
        }
        return counts;
    }
    
    public static int[][] CriticalITCount(Complex g2, Complex g3, int n,
        int bail, double tol, int maxOrbitLength) {
        
        Complex[] criticalValues = wpCriticalValues(g2,g3);
        int i,k;
        Complex[] recentOrbit;
        double[] currentDist;
        int[][] result = new int[3][2];
        
        for(k=0; k<3; k++) {
			i = 0;
			recentOrbit = new Complex[maxOrbitLength];
			Complex z = criticalValues[k];
			while(i < maxOrbitLength) {
				z = wp(z, g2, g3, n);
				recentOrbit[i] = z;
				i++;
			}
			
			z = wp(z, g2, g3, n);
			currentDist = dist(z, recentOrbit);
			recentOrbit[i % maxOrbitLength] = z;
			while(i<bail && currentDist[0]>tol) {
				z = wp(z, g2, g3, n);
				currentDist = dist(z, recentOrbit);
				i = i+1;
				recentOrbit[i % maxOrbitLength] = z;
			}
			result[k][0] = i;
			if((i % maxOrbitLength) > currentDist[1]) {
				result[k][1] = i % maxOrbitLength - (int) currentDist[1];
			}
			if((i % maxOrbitLength) <= currentDist[1]) {
				result[k][1] = maxOrbitLength - ((int) currentDist[1] - i % maxOrbitLength);
			}
        }
		return result;
    }            



    public static double[] CriticalITCountSquare(Complex g2, int n,
        int bail, double tol, int maxOrbitLength) {
    
        int i;
        Complex[] recentOrbit = new Complex[maxOrbitLength];
        double[] currentDist;
        Complex zeroComplex = new Complex(0.0,0.0);
        double maxSize = 0;
        double maxSizeTest;
        double[] result = new double[3];
        
        i = 0;
        Complex z = g2.sqrt().scale(0.5);
        while(i < maxOrbitLength) {
            z = wp(z, g2, zeroComplex, n);
            maxSizeTest = z.abs();
            if(maxSizeTest > maxSize && i < 15) {
                maxSize = maxSizeTest;
            }
            recentOrbit[i] = z;
            i++;
        }
        
        z = wp(z, g2, zeroComplex, n);
        currentDist = dist(z, recentOrbit);
        recentOrbit[i % maxOrbitLength] = z;
        int  out1 = i % maxOrbitLength;
        while(i<bail && currentDist[0]>tol) {
            z = wp(z, g2, zeroComplex, n);
            maxSizeTest = z.abs();
//             if(maxSizeTest > maxSize) {
//                 maxSize = maxSizeTest;
//             }
            currentDist = dist(z, recentOrbit);
            i = i+1;
            recentOrbit[i % maxOrbitLength] = z;
        }
        result[0] = (double) i;
        if((i % maxOrbitLength) > currentDist[1]) {
            result[1] = (double) i % maxOrbitLength - currentDist[1];
        }
        if((i % maxOrbitLength) <= currentDist[1]) {
            result[1] = (double) maxOrbitLength - (currentDist[1] - (double) i % maxOrbitLength);
        }
        result[2] = (double) maxSize;
        return result;
    }


    // To compute g2-space, we iterate P(z; g2, 0) starting from e1 = sqrt(g2)/4.
    // Note that e2 = -sqrt(g2)/4 has the same orbit as e1 and e3=0 is in the
    // the Julia set.  Thus, P(z;g2,0) can have at most one attractive orbit and
    // the question is whether e1 converges to such an orbit or not.
    public static double[][][] g2Space(Complex g2Min, Complex g2Max,
        int n, int bail, double tol, int resRe, int resIm, int maxOrbitLength) {
        
        Complex g2 = new Complex();
        int j, k;
        double[][][] counts = new double[resIm + 1][resRe + 1][3];
        
        for(k=0; k <= resIm; k++) {
            for(j=0; j <= resRe; j++) {
                g2 = new Complex((double) (g2Max.re() - g2Min.re())*j/resRe + g2Min.re(),
                    (double) (g2Max.im() - g2Min.im())*k/resIm + g2Min.im());
                counts[k][j] = CriticalITCountSquare(g2, n, bail, tol, maxOrbitLength);
            }
        }
        return counts;
    }
    
    public static double[][] g2Triangle(double radius, int n, int bail, double tol,
        int res, int maxOrbitLength) {
        int i, j, cnt;
        Complex g2;
        double[][] counts = new double[(res+1)*(res+2)/2][3];
        
        cnt = 0;
        for(i=0; i <= res; i++) {
            for(j=0; j <=res - i; j++) {
                g2 = new Complex((double) i*radius/res + j*radius/(2*res),
                    (double) j* 0.8660254037844386*radius/res);
                counts[cnt++] = CriticalITCountSquare(g2, n, bail, tol, maxOrbitLength);
            }
         }
         return counts;
     }

    public static int[][][][] g2ParallelSpace(Complex g2Min, Complex g2Max, Complex g3,
        int n, int bail, double tol, int resRe, int resIm) {
        
        Complex g2 = new Complex();
        int j, k;
        int[][][][] counts = new int[resIm + 1][resRe + 1][3][2];
        
        for(k=0; k <= resIm; k++) {
            for(j=0; j <= resRe; j++) {
                g2 = new Complex((double) (g2Max.re() - g2Min.re())*j/resRe + g2Min.re(),
                    (double) (g2Max.im() - g2Min.im())*k/resIm + g2Min.im());
                counts[k][j] = CriticalITCount(g2, g3, n, bail, tol, 10);
            }
        }
        return counts;
    }

}