Complexity of the Mandelbrot set on rationals

Question

Given two rationals $a,b \in \mathbb{Q}$, call $c = a + ib$, i.e., the complex number represented by these two rationals.

A point $c$ is contained within the Mandelbrot set $M$ if the following procedure never halts:

$z = c$

$while (abs(z) <= 2)$

$\hspace{0.3in} z = z^2+c$

Normally we pick some $k$ (say, 50) and then if it doesn't halt after that many iterations we assume it is in the Mandelbrot set and stop looping. The ordinary algorithm for this takes $O(n2^k)$ time (simply computing the above loop using something like Java's BigInteger class for the numerators and denominators and using fraction arithmetic).

Is there a polynomial time algorithm for determining if a given $c$ breaks out of this loop within $k$ steps, in terms of the magnitude of $k$ and in terms of $n$ bits representing the numerator and denominator of $a$ and $b$?

For reference, $abs(c) = \sqrt{a^2 + b^2}$ and $c^2 = a^2 + 2abi +b^2i^2 = a^2 - b^2 + 2abi$ because $i^2 = -1$.

(this was cross posted at cstheory stackexchange but in hindsight I thought it actually fit better here)

Answer 1

You probably cannot practically do 100 iterations, at least not in the current universe. Let $$ f_c(z) = z^2 + c $$ and write $f_c^{\circ n}(z)$ for the $n$'th iterate of $f_c$. You want to start with a Gaussian rational number $\alpha\in\mathbb Q[i]$ and compute the exact value of $f_c^{\circ n}(\alpha)$. (I realize that you're only interested in taking $\alpha=c$, but the complexity doesn't really change that much.) Let's write the number $f_c^{\circ n}(\alpha)$ as $$f_c^{\circ n}(\alpha) = \frac{A_n}{B_n}+\frac{C_n}{D_n}i,$$ where $A_n,B_n,C_n,D_n\in\mathbb Z$ and the fractions are in lowest terms. The problem is that the number of digits in these integers grows enormously rapidly. To be precise, the following limit converges: $$ \hat h_c(\alpha) := \lim_{n\to\infty} \frac{\log\max\{|A_n|,|B_n|,|C_n|,|D_n|\}}{2^n}, $$ and it has the property that $$ \hat h_c(\alpha) > 0 \quad\text{if and only if $\alpha$ is not perperiodic.}$$ So if your want to compute $f_c^{\circ n}(\alpha)$ exactly for a non-perperiodic $\alpha$, you'll need to work with integers having roughly $2^n$ digits! If $n=100$, this would seem to be infeasible.

If you want more information about $ \hat h_c(\alpha)$, you can google dynamical canonical height.

Answer 2

Okay so I haven't proved this yet but in practice this seems to be working pretty well for me. The trick is that you use interval arithmetic for the real and imag component of your number, using fractions with biginteger type numbers for their numerator and denominator. The bits still grow very quickly, but what you can do is "bound" them at each iter. For example, 0.000123141532 can become [0.000124, 0.000123], then you use those as your new intervals. If you end up having your upper estimate of the norm get >= 4 but the lower estimate of the norm be < 4 then you want to try a closer interval (so like [0.00012314, 0.00012315]). I find that usually 200 bits is plenty enough precision to run 1000 iters which wouldn't be practical otherwise, but these could grow as much as needed.

Of course this code isn't super fast and optimized and using floating point numbers for rendering fractals is much easier, but the important part here is that these results are exact, this is guaranteed to give you whether or not it breaks in k steps. It runs basically in time linear in k since you can just fix precision, but there may be worst case behaviors on the edge of the mandelbrot that require the interval to get so small that you get exponential slowdown. There was some stuff on perturbation theory that might help prove that this doesn't happen but I haven't done that yet. I haven't ran into that issue in all my renders I've done with this so far at least.

For now, here is that idea implemented in python. I'd like to clean up this code/comment it but I figured I'd just share it first so others can believe me that this works and maybe help with proving perturbation theory stuff I'm not sure about.

Specifically, here is the interval class I made in python:

import math
import fractions
class Interval(object):
  def __init__(self, lo, hi=None):
    if type(lo) is Interval and hi is None:
      self.lo = lo.lo
      self.hi = lo.hi
    else:
      self.lo = fractions.Fraction(lo)
      if hi is None:
        hi = lo
      self.hi = fractions.Fraction(hi)
    if self.lo > self.hi:
      raise Exception("lo " + str(self.lo) + " is greater than hi: " + str(self.hi))

  def applyOp(self, val, op):
    val = Interval(val)
    loop = getattr(self.lo, op)
    hiop = getattr(self.hi, op)
    things = [loop(val.lo), loop(val.hi), hiop(val.lo), hiop(val.hi)]
    return Interval(min(things), max(things))

  def __mul__(self, val):
    return self.applyOp(val, '__mul__')

  def __rmul__(self, val):
    return self.applyOp(val, '__rmul__')

  def __add__(self, val):
    return self.applyOp(val, '__add__')

  def __sub__(self, val):
    return self.applyOp(val, '__sub__')

  def __div__(self, val):
    return self.applyOp(val, '__div__')

  def __neg__(self):
    return Interval(-self.hi, -self.lo)

  def __repr__(self):
    return self.__str__()

  def __pow__(self, val):
    if val == 2:
      if self.lo <= 0 and self.hi >= 0:
        return Interval(0, max(self.lo*self.lo, self.hi*self.hi))
      elif self.lo <= 0 and self.hi <= 0:
        return Interval(self.hi*self.hi, self.lo*self.lo)
      elif self.lo >= 0 and self.hi >= 0:
        return Interval(self.lo*self.lo, self.hi*self.hi)

    raise Exception("don't know how to do higher pow yet")

  def reduce(self, maxNDigits):
    lower = reduceFract(self.lo, maxNDigits, roundDown=True)
    bigger = reduceFract(self.hi, maxNDigits, roundDown=False)
    return Interval(lower, bigger)

  def __str__(self):
    res = "[" + str(self.lo) + "," + str(self.hi) + ']'
    return res

  def approx(self):
    return [float(self.lo), float(self.hi)]

  def __contains__(self, other):
    other = Interval(other)
    return self.lo <= other.lo and self.hi >= other.hi

def reduceFract(res, biggestAllowed, roundDown=True):
  if res.numerator == 0 or res.denominator == 0: return res
  mag = int(min(math.log(abs(res.numerator), 2),math.log(abs(res.denominator), 2)))
  if mag > biggestAllowed:
    mag = mag - biggestAllowed
    num = res.numerator>>mag
    den = res.denominator>>mag
    resFrac = fractions.Fraction(num, den)
    if roundDown and resFrac > res:
      resFrac = fractions.Fraction(num-1, den)
    elif not roundDown and resFrac < res:
      resFrac = fractions.Fraction(num+1, den)
    res = resFrac
  return res

now we can do this:

def mandelInterval(cr, ci, n, reduceAmount):
  cr = Interval(cr, cr)
  ci = Interval(ci, ci)
  resr, resi, state = mandelIntervalHelper(cr, ci, n, reduceAmount)
  mag = resr*resr + resi*resi
  return resr.approx(), resi.approx(), mag.approx(), state

def mandelIntervalHelper(cr, ci, n, reduceAmount):
  if n <= 0:
    return cr, ci, None
  else:
    i = 1
    pr, pi, state = mandelIntervalHelper(cr,ci,n-1, reduceAmount)
    while True:
      if state == "need better estimates":
        i += 1
        pr, pi, state = mandelIntervalHelper(cr,ci,n-1, reduceAmount*i)
      else:
        break
    if state in ['outside', 'need better estimates']: return pr, pi, state
    resr = (pr*pr-pi*pi + cr).reduce(reduceAmount)
    resi = (2*pr*pi+ci).reduce(reduceAmount)
    mag = resr**2 + resi**2
    if mag.lo >= 4 and mag.hi >= 4:
      return resr, resi, "outside"
    elif mag.lo < 4 and mag.hi >= 4:
      return resr, resi, "need better estimates"
    return resr, resi, "inside"     

So for example, saving this to ms.py:

>>> import ms
>>> # c= -0.397959183673 + -0.673469387755i
>>> # 1000 iters, 20 bits of precision (goes to 40, 60, 80, ... if needed)
>>> ms.mandelInterval(-0.397959183673, -0.673469387755, 1000, 20)
([-0.05965659022331238, -0.033619076013565063], [-2.383112907409668, -2.357083320617676], [5.556972022606055, 5.682786038219633], 'outside')
>>> # [-0.05965659022331238, -0.033619076013565063] is interval for real(c)
>>> # [-2.383112907409668, -2.357083320617676] is interval for imag(c)
>>> # [5.556972022606055, 5.682786038219633] is interval for distance
>>> # because both sides of the interval for distance are >= 4, the last value is 'outside', meaning we are outside the mandelbrot set. This will be 'inside' if we are inside it.