import java.util.*;
import java.awt.image.*;

public class MandelbrotSet 
{
	// This variant requires a power of 2
	public int w = 512;
	public int h = 512;
	public double xs = -2.0;
	public double ys = -2.0;
	public double d = 4.0;
	public int stepsMax = 200;
	int colors[];
	double[] zp;
	double[] zip;
//	DoubleDouble[] ddzp;
//	DoubleDouble[] ddzip;
	public boolean computing = false;
	public boolean stopRequest = false;
	public int pix[];
	int iterations[];
	int[] rx1;
	int[] ry1;
	int[] rx2;
	int[] ry2;
	int[] px;
	int[] py;
	MandelbrotUser user;
	MandelbrotSet(MandelbrotUser userArg,
		int wArg, int hArg, double xsArg, double ysArg, 
		double dArg, int stepsMaxArg)
	{
//		System.out.println("w: " + wArg + " h: " + hArg + " xs: " +
//			xsArg + " ys: " + ysArg + " d: " + dArg + " i: " +
//			stepsMaxArg);
		user = userArg;
		w = wArg;
		h = hArg;
		pix = new int[w * h];
		iterations = new int[w * h];
		rx1 = new int[w * h];
		ry1 = new int[w * h];
		rx2 = new int[w * h];
		ry2 = new int[w * h];
		px = new int[16 + w * 4];
		py = new int[16 + w * 4];
		xs = xsArg;
		ys = ysArg;
		d = dArg;
		stepsMax = stepsMaxArg;
		newSteps();
	}
	int orbitsFound = 0;
	void renderSet()
	{
		computing = true;
		int x, y;
		int low = -1, high = -1;
		int ptotal = w * h;
		for (x = 0; (x < ptotal); x++) {
			pix[x] = 0xFF000000;
			iterations[x] = -1;
		}	 	
		int inc = w;
		int rcount = 0;
		int[] rx = rx1;
		int[] ry = ry1;
		int[] rxn = null;
		int[] ryn = null;
		rx[rcount] = 0;
		ry[rcount] = 0;
		rcount++;
		int shift = 2;
		int tsubs = 16;
		int pass = 0;
		int calculated = 0;
		int orbitsFound = 0;
//		long start = System.currentTimeMillis();
		while (inc > 1) {
			// Allow only 4 points as a hack to allow
			// resolutions that are powers of 2 but not of 4
			if (inc == 2) {
				shift = 1;
				tsubs = 4;
			}
			int r;
			if (rx == rx1) {
				rxn = rx2;
				ryn = ry2;
			} else {
				rxn = rx1;
				ryn = ry1;
			}
			int rcountn = 0;
			for (r = 0; (r < rcount); r++) {
				int lx = rx[r];		
				int ly = ry[r];		
				int xi, yi;
				int n = 0;
				for (yi = ly; (yi < (ly + inc)); yi += (inc >> shift)) {
					for (xi = lx; (xi < (lx + inc)); xi += (inc >> shift)) {
						px[n] = xi;
						py[n] = yi;
						n++;
					}
				}	
				for (yi = ly; (yi < (ly + inc)); yi += 2) {
					px[n] = lx;
					py[n] = yi;
					n++;
					px[n] = lx + inc - 1;
					py[n] = yi;
					n++;
				}
				for (xi = lx; (xi < (lx + inc)); xi += 2) {
					px[n] = xi;
					py[n] = ly;
					n++;
					px[n] = xi;
					py[n] = ly + inc - 1;
					n++;
				}
				int p;
				for (p = 0; (p < n); p++) {
					x = px[p];	
					y = py[p];	
					if (iterations[y * w + x] != -1) {
						continue;
					}
					int i;
					int steps = -1;
					// Iteration is vastly more expensive
					// than this little boolean test, we
					// can afford to do it more often
					// for better interactive response
					if (stopRequest) {
						computing = false;
						return;
					}	
					steps = doubleIterate(x, y);
//					steps = doubleDoubleIterate(x, y);
					if (steps != -1) {
						iterations[y * w + x] = steps;
						if ((low == -1) || (steps < low)) {
							low = steps;
						}
						if ((high == -1) || (steps > high)) {
							high = steps;
						}
					} else {
						// Don't redo the set itself
						iterations[y * w + x] = -2;
					}
					calculated++;
				}
				int it = iterations[py[0] * w + px[0]];
				for (p = 1; (p < n); p++) {
					if (iterations[py[p] * w + 
						px[p]] != it) {	
						break;
					}
				}
				if (p != n) {
					// tsubs is right here
					for (p = 0; (p < tsubs); p++) {
						rxn[rcountn] = px[p];	
						ryn[rcountn] = py[p];	
						rcountn++;
					}
				} else {
					int sx, sy;
					for (sy = ly; (sy < ly + inc); sy++) {	
						for (sx = lx; (sx < lx + inc); sx++) {	
							iterations[sy * w + sx] = iterations[ly * w + lx];
						}
					}
				}
			}
			int i;
			if ((low != -1) && (high != -1)) {
				for (i = low; (i < (high + 1)); i++) {
					int is;
					int range = (high + 1 - low) / 5;
					if (range == 0) {
						range = 1;
					} 
					int segment = ((i - low) / range);
					int offset = ((i - low) % range);
					int redi, greeni, bluei;
					switch(segment) {
						case 0:
						redi = 255 - offset * 255 / range;	
						greeni = 0;
						bluei = 255;
						break;
						case 1:
						redi = 0;
						greeni = offset * 255 / range;
						bluei = 255;
						break;
						case 2:
						redi = 0;
						greeni = 255;
						bluei = 255 - offset * 255 / range;
						break;
						case 3: 
						redi = offset * 255 / range;
						greeni = 255;
						bluei = 0;
						break;
						case 4:
						redi = 255;
						greeni = 255 - offset * 255 / range;
						bluei = 0;
						break;
						default:
						redi = 255;
						greeni = 0;
						bluei = 0;
						break;
					}
		//			System.out.println(redi + "," + greeni + "," + bluei);
					colors[i] = 0xFF000000 + (redi << 16) + (greeni << 8) + bluei;
				} 	
				int sub = inc >> shift;
				for (y = 0; (y < h); y += sub) {
					if (stopRequest) {
						computing = false;
						return;
					}	
					for (x = 0; (x < w); x += sub) {
						int it = iterations[y * w + x];
						int color = 0xFF000000;	
						if (it >= 0) {
							color = colors[it];
						}
						int xi, yi;
						for (yi = y; ((yi < y + sub) && (yi < h)); yi++) {
								for (xi = x; ((xi < x + sub) && (xi < w)); xi++) {
										pix[yi * w + xi] = color;
								}
						}
					}
				}
			}
			// TODO PICKUP HERE
			if (pass > 1) {
				if (user != null) {
					user.MandelbrotUpdated();
				}
			}	
			inc >>= shift;
			rx = rxn;
			ry = ryn;
			rcount = rcountn;
			pass++;
		}
//		long end = System.currentTimeMillis();
//		System.out.println("elapsed ms: " + (end - start));
//		System.out.println(calculated + " of " + ptotal + ": " +
//			(double) calculated * 100.0 / (double) ptotal + ", orbits: " +
//			orbitsFound + ", " +
//			(double) orbitsFound * 100.0 / (double) ptotal);
		computing = false;
		if (user != null) {
			user.MandelbrotFinished();
		}
	}
	void newSteps()
	{
		colors = new int[stepsMax];
		zp = new double[stepsMax];
		zip = new double[stepsMax];
//		ddzp = new DoubleDouble[stepsMax];
//		ddzip = new DoubleDouble[stepsMax];
//		int i;
//		for (i = 0; (i < stepsMax); i++) {
//			ddzp[i] = new DoubleDouble(0);
//			ddzip[i] = new DoubleDouble(0);
//		}
	}
	int doubleIterate(int x, int y)
	{
		int steps = -1;
		double c = (((double) x / (double) w) * d) + xs;
		double ci = (((double) y / (double) w) * d) + ys;
		// Use orbit detection
		double zf = 0;
		double zfi = 0;
		double zs = 0;
		double zsi = 0;
		int i;
		for (i = 0; (i < stepsMax); i += 2) {
			zp[i] = zf;
			zip[i] = zfi;
			double tzi = (zf * zfi) * 2;
			zf = (zf * zf) - (zfi * zfi);
			zfi = tzi;
			zf += c;	
			zfi += ci;	
			zp[i] = zf;
			zip[i] = zfi;
	
			if (zf * zf + zfi * zfi > 4.0) {
				steps = i;
				break;
			}		 	 
			
			tzi = (zf * zfi) * 2;
			zf = (zf * zf) - (zfi * zfi);
			zfi = tzi;
			zf += c;	
			zfi += ci;	
	//						tzi = (zs * zsi) * 2;
	//						zs = (zs * zs) - (zsi * zsi);
	//						zsi = tzi;
	//						zs += c;	
	//						zsi += ci;	
			if (zf * zf + zfi * zfi > 4.0) {
				steps = i + 1;
				break;
			}		 	 
			// Orbit detected
			if ((i > 0) && (zf == zp[i >> 1]) && (zfi == zip[i >> 1])) {
				orbitsFound++;
				steps = -1;
				break;
			}
		}
		return steps;
	}
//	// First pass. Still todo: the parameters should also be DoubleDouble,
//	// more calculations should be DoubleDouble, more variables should
//	// persist, text output of DoubleDouble still broken.
//	int doubleDoubleIterate(int x, int y)
//	{
//		int steps = -1;
//		DoubleDouble ddx = new DoubleDouble(x);
//		DoubleDouble ddy = new DoubleDouble(y);
//		DoubleDouble ddxs = new DoubleDouble(xs);
//		DoubleDouble ddys = new DoubleDouble(ys);
//		DoubleDouble ddd = new DoubleDouble(d);
//		DoubleDouble quot = new DoubleDouble();
//		DoubleDouble ddw = new DoubleDouble(w);
//		quot.div(ddx, ddw);
//		quot.mult(quot, ddd);
//		DoubleDouble c = new DoubleDouble();
//		c.add(quot, ddxs);
//		quot.div(ddy, ddw);
//		quot.mult(quot, ddd);
//		DoubleDouble ci = new DoubleDouble();
//		ci.add(quot, ddys);
////		System.out.println(c + "," + ci);	
//		// Use orbit detection
//		DoubleDouble zf = new DoubleDouble(0);
//		DoubleDouble zfi = new DoubleDouble(0);
//		DoubleDouble zs = new DoubleDouble(0);
//		DoubleDouble zsi = new DoubleDouble(0);
//		DoubleDouble tzi = new DoubleDouble(0);
//		DoubleDouble prod = new DoubleDouble(0);
//		DoubleDouble prod2 = new DoubleDouble(0);
//		int i;
////		System.out.println("beginning steps");
//		for (i = 0; (i < stepsMax); i += 2) {
////			System.out.println(zf + "," + zfi);
//			ddzp[i].set(zf);
//			ddzip[i].set(zfi);
//			tzi.mult(zf, zfi);
//			tzi.mult(tzi, 2.0);
//			zf.mult(zf, zf); 
//			prod.mult(zfi, zfi);
//			prod.negate();
//			zf.add(zf, prod);
//			zfi.set(tzi);
//			zf.add(zf, c);
//			zfi.add(zfi, ci);
////			System.out.println("zf and zfi step " + i + ": " +
////				zf.hi + "," + zf.lo + " " + zfi.hi + "," +
////				zfi.lo); 
//			ddzp[i].set(zf);
//			ddzip[i].set(zfi);
//			prod.mult(zf, zf);	
//			prod2.mult(zfi, zfi);
//			prod.add(prod, prod2);
//			if (prod.gt(prod, 4.0)) {
////				System.out.println("product too large: " + prod);
//				steps = i;
//				break;
//			}
//			tzi.mult(zf, zfi);
//			tzi.mult(tzi, 2.0);
//			prod.mult(zf, zf);
//			prod2.mult(zfi, zfi);
//			prod2.negate();
//			zf.add(prod, prod2);
//			zfi.set(tzi);
//			zf.add(zf, c);
//			zfi.add(zfi, ci);
//			prod.mult(zf, zf);
//			prod2.mult(zfi, zfi);
//			prod.add(prod, prod2);
//			if (prod.gt(prod, 4.0)) {
//				steps = i + 1;
////				System.out.println("product too large: " + prod);
//				break;
//			}		 	 
//			// Orbit detected
//			if ((i > 0) && zf.eq(ddzp[i >> 1]) && zfi.eq(ddzip[i >> 1])) {
//				orbitsFound++;
//				steps = -1;
//				break;
//			}
//		}
////		System.out.println(steps);
//		return steps;
//	}
}