Buddhabrot Buddhabrot
Interior and exterior distance bounds for the Mandelbrot

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?

Interior distance bounding

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);
  }
Exterior distance bounding

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.

 
Copyright © Albert Lobo Cusidó 2006-2014