My BBrot zooming technique involves estimating the distance from a point C to the nearest point on the boundary of the Mandelbrot set. Since these distance bounding methods are fairly complicated and non-intuitive, I found it appropriate to dedicate this appendix to them. Also, very few explanations are available on the internet (especially for the interior bound); as expected, MuEncy comments on both bounding methods, and even gives souce for the exterior bounder.
Offline, you can look up a book named "The science of fractal images", by Heinz-Otto Peitgen & Dietmar Saupe editors. This book is worth its weigh in gold -out-of-print, I reckon-. Why doesn't anyone upload it as an e-Book?
The first thing you have to remember is that the resulting distances should not be considered precise. They are only bounds; the closest point on the M border is never found. Well, that's because of the nature of the algorithm.
What we'll be doing, plain and simple, is looking for disks completely contained in M. Maybe a bigger disk would still be inside M, but that doesn't matter, because we're just calculating bounds. The algorithm relies on a theorem by Hubbard & Douady (who demonstrated that M is connected by mapping M onto the unit disk), combined with the Koebe 1/4 theorem.
The process, if I gathered correctly, is as follows:
- Map C onto the unit disk using Hubbard & Douady's method.
- Get the inverse of the map using the Möbius transformation.
- Using Hubbard & Douady's method again, we "unmap" from the inverse map to get a complex value. Its magnitude, d, will be used to compute a lower bound.
- The Koebe 1/4 theorem applies, by construction; hence, all points within a radius of d/4 from C are sure to be in M.
In practice, not only do we need C to be inside M, but we also require C's periodic cycle. That is increasingly harder to detect as C approaches the M border. The following Java snippet is an example of how to find the cycle.
/** * Finds the periodic cycle for a given C, provided C is in M, * and the number of iterates to check is sufficient. * @param Cr double - Real component of C. * @param Ci double - Imaginary component of C. * @param maxN int - Maximum number of iterates to check. * @return double[][] - The periodic cycle. The period equals * 'getPeriodicCycle(Cr,Ci,maxN).length'. Returns null if C * is outside M, and an empty array if the periodic cycle was not found. */ public double[][] getPeriodicCycle(double Cr, double Ci, int maxN) { double Zr, Zi, ZrT, c1, c2; double[] iterR = new double[maxN]; double[] iterI = new double[maxN]; int i, j, n, maxNT; // We'll do various checks to detect more obvious cycles faster. int[] nsToChk = new int[] { 0, 0, 0, 0, maxN }; for (i = 3; i > -1; i--) nsToChk[i] = nsToChk[i + 1] / 2; // Start iterating. for (n = 1, iterR[0] = Zr = iterI[0] = Zi = i = 0; i < 5; i++) { maxNT = nsToChk[i]; while (n < maxNT) { iterR[n] = ZrT = Zr * Zr - Zi * Zi + Cr; iterI[n++] = Zi = 2.0 * Zr * Zi + Ci; Zr = ZrT; // If C is outside M, return null. if (Zr * Zr + Zi * Zi > 4) return null; } // This is the periodicity check. for (j = n - 2; j > -1; j--) { c1 = Zr - iterR[j]; c2 = Zi - iterI[j]; if (c1 * c1 + c2 * c2 < 1E-30) { // Okay, close enough. int period = n - j - 1; double[][] cycle = new double[period][2]; int offset = ((n / period) - 1) * period; for (int k = 0; k < period; k++, offset++) { cycle[k][0] = iterR[offset]; cycle[k][1] = iterI[offset]; } return cycle; } } } // Iterates didn't escape, but period never became apparent // -return empty array. return new double[0][0]; }
Once we have the periodic cycle, the distance estimation is pretty straight-forward -see the Java source below:
/** * Calculates a lower bound for the distance between C * and the closest point in the border of M. * @param cycle double[][] - C's periodic cycle. * @return double - The interior distance estimation. */ public double getInteriorLowerBound(double[][] cycle) { // Real and imaginary components for complex numbers D1, D2, D3, and D4; // and temporary variables to store the complex numbers in recursion. double Zr, Zi, D1r, D1i, D2r, D2i, D3r, D3i, D4r, D4i; double D1rT, D1iT, D2rT, D2iT, D3rT, D3iT, D4rT, D4iT; int i, period = cycle.length; // Initial values: D1 = 1; D2 = 0; D3 = 0; D4 = 0; D1r = 1; D1i = D2r = D2i = D3r = D3i = D4r = D4i = 0; // Start iterating. for (i = 0; i < period; i++) { // No need to iterate Z; values are already in 'cycle'. Zr = cycle[i][0]; Zi = cycle[i][1]; // D1 = 2 * Z * D1; D1rT = 2 * (Zr * D1r - Zi * D1i); D1iT = 2 * (Zi * D1r + Zr * D1i); // D2 = 2 * Z * D2 + 1; D2rT = 2 * (Zr * D2r - Zi * D2i) + 1; D2iT = 2 * (Zi * D2r + Zr * D2i); // D3 = 2 * (D1^2 + Z * D3); D3rT = 2 * ((Zr * D3r - Zi * D3i) + (D1r * D1r - D1i * D1i)); D3iT = 2 * ((Zr * D3i + Zi * D3r) + (2 * D1r * D1i)); // D4 = 2 * (D1 * D2 + Z * D4); D4rT = 2 * ((Zr * D4r - Zi * D4i) + (D1r * D2r - D1i * D2i)); D4iT = 2 * ((Zr * D4i + Zi * D4r) + (D1r * D2i + D1i * D2r)); // Update variables. D1r = D1rT; D1i = D1iT; D2r = D2rT; D2i = D2iT; D3r = D3rT; D3i = D3iT; D4r = D4rT; D4i = D4iT; } // A = 1 - |D1|^2; double A = 1 - (D1r * D1r + D1i * D1i); // B = |D4 + D3 * (D2 / (1 - D1))|; double B = (1 - D1r) * (1 - D1r) + D1i * D1i; D1rT = (D2r * (1 - D1r) - D2i * D1i) / B; D1iT = (D2i * (1 - D1r) + D2r * D1i) / B; D2rT = D4r + (D3r * D1rT - D3i * D1iT); D2iT = D4i + (D3i * D1rT + D3r * D1iT); B = StrictMath.sqrt(D2rT * D2rT + D2iT * D2iT); // Return lower bound -that is, 1/4 the estimated bound. return A / (4.0 * B); }
Computing a lower bound for the distance between a point C outside M and M is even more simple. As in the interior bounding, we rely on theorems by Hubbard & Douady, and Koebe. The actual algorithm is based on the ideas of J. Milnor and W. Thurston, expanded by A. Wilkes & J. Hubbard. Beware, though, the algorithm works well for Cs reasonably close to M.
This is the Java source:
/** * Calculates a lower bound for the distance between C and the closest * point in the border of M, if C is outside M. * @param Cr double - Real component of C. * @param Ci double - Imaginary component of C. * @param maxN int - Maximum number of times the Mandelbrot function will be iterated. * @return double - The approximate distance between C and the closest * point in M. If C is in M, -1 will be returned. */> public double getExteriorLowerBound(double Cr, double Ci, int maxN) { double Zr, Zi, ZrT, ZiT, c1, c2; double[] iterR = new double[maxN]; double[] iterI = new double[maxN]; int i, j, n, maxNT; // We'll do various checks to detect periods faster. int[] nsToChk = new int[] { 0, 0, 0, 0, maxN }; for (i = 3; i > -1; i--) nsToChk[i] = nsToChk[i + 1] / 2; // Start iterating. for (n = 1, iterR[0] = Zr = iterI[0] = Zi = i = 0; i < 5; i++) { maxNT = nsToChk[i]; while (n < maxNT) { iterR[n] = ZrT = Zr * Zr - Zi * Zi + Cr; iterI[n++] = Zi = 2.0 * Zr * Zi + Ci; Zr = ZrT; // If C is outside M, do distance estimation. if (Zr * Zr + Zi * Zi > 4) { double dZr, dZi, dZrT; int nT = n - 1; // Recursively calculate first derivative: D = 2 * Z * D + 1; for (dZr = dZi = i = 0; i < nT; i++) { ZrT = iterR[i]; ZiT = iterI[i]; dZrT = 2.0 * (ZrT * dZr - ZiT * dZi) + 1; dZi = 2.0 * (ZrT * dZi + ZiT * dZr); dZr = dZrT; } // Lower bound: (2 * log |Z| * (|Z| / |D|)) / 4; c1 = StrictMath.sqrt(Zr * Zr + Zi * Zi); c2 = StrictMath.sqrt(dZr * dZr + dZi * dZi); return (StrictMath.log(Zr * Zr + Zi * Zi) * c1 / c2) / 4.0; } } // Periodicity check. for (j = n - 2; j > -1; j--) { c1 = Zr - iterR[j]; c2 = Zi - iterI[j]; if (c1 * c1 + c2 * c2 < 1E-30) return -1.0; } } return -1.0; }
Both interior and exterior bounding can be combined to produce images such as this:
For lack of a better color palette, I used gray-scale.
Points outside M get brighter as they approach the border. Points in M get darker towards the border.