# Fast External Ray Algorithms
There are two related algorithms that relate the symbolic properties of the Mandelbrot set (patterns of periods, angled internal addresses and external angles) with its numeric properties (coordinates in the complex plane, i.e. \(c \in \mathbb{C}\)):
- tracing external rays inwards to find \(c\) from \(\theta\)
- tracing external rays outwards to find \(\theta\) from \(c\)
# 1 Problem
However simple implementation is prohibitively expensive: to trace a ray to depth \(O(p)\) requires \(O(p^2)\) multiplications (because you need \(O(p)\) steps, with linearly increasing cost) each of which costs \(M(p)\) (the cost of multiplying numbers with precision proportional to \(p\)). The simplest multiplication algorithm has \(M(p) = O(p^2)\), but libraries using Karatsuba algorithm and so on get below \(M(p) = O(p^{1.58})\). Thus the simplest implementation is \(O(p^4)\), with advanced multiplication \(O(p^{3.58})\).
# 2 Solution
Good news: for some external angles, like those that arise in practice via artistic zooming methods like Julia morphing, this can be improved using perturbation techniques. Finding a nucleus of period \(p\) given a good initial guess for Newton’s method costs \(O(p M(p)) \approx O(p^{2.58})\). Finding its orbit and rounding to fixed low precision costs the same. Then the cost of evaluating \(f_c^n(0)\) is reduced from \(O(p M(p))\) using full precision, to \(O(p)\) using low fixed precision with perturbation. Thus the cost of tracing the external ray is reduced from \(O(p^2 M(p))\) using full precision, to \(O(p^2)\) using low fixed precision with perturbation. The total cost is therefore \(O(p M(p))\), that is, it is dominated by the reference because \(M(p) > O(p)\). That is, the asymptotic cost is reduced by a factor of \(p\), assuming we have a nearby \(c\) to use as a reference!
# 3 Example
For example, consider JWM’s Old Wood Dish. It has angled internal address \[1 \frac12 2 \frac{16}{17} 33 \frac12 34 \frac13 69 \frac12 70 \cdots \frac13 (p-1) \frac12 p \frac13 (2p+1) \frac12 (2p+2) \cdots\] Each truncation of the address corresponds to a minibrot island, and the next in the sequence is relatively nearby. The external angles follow a regular pattern too: if the angles are denoted \(.(l_p)\) and \(.(h_p)\) where \(l\), \(h\) are blocks of binary digits, then \(l_p = l_{p-1}1\) and \(h_p = h_{p-1}0\), and \(l_{2p+1} = l_pl_p1\) and \(h_{2p+1} = l_ph_p0\).
# 4 Algorithm
Using Old Wood Dish example above; other locations with different patterns in the angled internal addresses will need corresponding adjustments.
-
input \(c_{p-1}\) is a nucleus of period \(p-1\) with precision \(k(p-1)\).
-
trace external ray \(l_p\) to depth \(3p\) using \(c_{p-1}\) as reference.
-
use Newton’s method to find \(c_p\) from the endpoint of the ray with precision \(k(p)\).
-
trace external ray \(l_{2p+1}\) to depth \(3(2p+1)\) using \(c_p\) as reference.
-
use Newton’s method to find \(c_{2p+1}\) from the endpoint of the ray with precison \(k(2p+1)\).
-
set \(p-1 := 2p+1\) and repeat from the start.
The function \(k(p)\) to give the required precision from the period may depend on the location, though I suspect it to be \(O(p)\) in general.
The total asymptotic cost of this algorithm is \(O(p M(k(p)))\) which is roughly \(O(p^{2.58})\) assuming linear precison increase and Karatsuba multiplication.
# 5 Results
The practical cost drops from almost 5 hours, to 22 minutes to find the period 9214 nucleus in the series, with 8 of those minutes for the step from 9213 to 9214. I used 1024 bits of precision for all the high precision calculations for simplicity. The final coordinates found (with the radius being 10 times the atom domain size to the power of 1.125) are:
re = "-7.621179376797638683152412470302243466863693891184035754206961666662148701906112386691681679369649115315407442420538767146892200125937933980945482470837951427211150657111520062340131401482067699893352032301169768087501364329395853845097310808420664970247127366095008899576141493001623616618339359179484065675671e-01"
im = "-9.573175593807915218110127195889290086484763586266554107386326351819880557832346419854569617721368084312027345073397513423650480579076128444209078562800246183549293915393520482765375446620329755078370285429051961341356129731430766990443543286328234828509304655693149935499090847158289064145468761037805221829293e-02"
radius = "3.24832872644961e-188"