>>> 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 'done''outside', meaning we are outside the mandelbrot set. This will be 'inside' if we are inside it.
>>> 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 'done', meaning we are outside the mandelbrot set. This will be 'inside' if we are inside it.
>>> 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.
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
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 'done', meaning we are outside the mandelbrot set. This will be 'inside' if we are inside it.
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
So for example:
# 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 'done', meaning we are outside the mandelbrot set. This will be 'inside' if we are inside it.
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
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 'done', meaning we are outside the mandelbrot set. This will be 'inside' if we are inside it.
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 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:
# 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 'done', meaning we are outside the mandelbrot set. This will be 'inside' if we are inside it.