# Complexity of the Mandelbrot set on rationals

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)

• Pickiing 50 is so 1980's. These days we pick 65535. Oct 12 '17 at 7:01
• @AndrejBauer lol true, though that is with floating points (which are somewhat accurate even zoomed in a lot without needing too much precision via using perturbation theory) but if you are forced to use rationals and not accept rounding error is there an easy way to do even like 100 iterations? Oct 12 '17 at 7:06
• I don't have an answer, but the following might help. It is not known whether the following problem is algorithmically decidable: Given $p, q \in \mathbb{Q}$ and $k \in \mathbb{N}$, determine whether the open circle with center at $p + q i$ and radius $2^{-k}$ intersects the Mandelbrot set. Oct 12 '17 at 9:22
• @AndrejBauer Do you have a reference for this? Oct 12 '17 at 21:22
• @FedericoPoloni: Hertling, P. (2005), Is the Mandelbrot set computable?. Mathematical Logic Quarterly, 51: 5–18. doi:10.1002/malq.200310124 Oct 13 '17 at 6:54

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.

• True, but might you be able to use the fact that all you need to know is if the norm squared of z is less than or equal to 4 to get around having to actually represent all 2^k bits? (k is the number of iterations and n is the number of bits in the rationals in the problem as posed but I see what you mean) Oct 12 '17 at 14:34
• @Silverman perperiodic? I guess you mean preperiodic? Nov 15 '17 at 20:29
• @coudy Oops, typo! Somehow it's easier to type perper than it is to type preper. Nov 15 '17 at 22:31

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 __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.