# 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. cc \in \mathbb{C}):

# 1 Problem

However simple implementation is prohibitively expensive: to trace a ray to depth O(p)O(p) requires O(p2)O(p^2) multiplications (because you need O(p)O(p) steps, with linearly increasing cost) each of which costs M(p)M(p) (the cost of multiplying numbers with precision proportional to pp). The simplest multiplication algorithm has M(p)=O(p2)M(p) = O(p^2), but libraries using Karatsuba algorithm and so on get below M(p)=O(p1.58)M(p) = O(p^{1.58}). Thus the simplest implementation is O(p4)O(p^4), with advanced multiplication O(p3.58)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 pp given a good initial guess for Newton’s method costs O(pM(p))O(p2.58)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 fcn(0)f_c^n(0) is reduced from O(pM(p))O(p M(p)) using full precision, to O(p)O(p) using low fixed precision with perturbation. Thus the cost of tracing the external ray is reduced from O(p2M(p))O(p^2 M(p)) using full precision, to O(p2)O(p^2) using low fixed precision with perturbation. The total cost is therefore O(pM(p))O(p M(p)), that is, it is dominated by the reference because M(p)>O(p)M(p) > O(p). That is, the asymptotic cost is reduced by a factor of pp, assuming we have a nearby cc to use as a reference!

# 3 Example

For example, consider JWM’s Old Wood Dish. It has angled internal address 112216173312341369127013(p1)12p13(2p+1)12(2p+2)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 .(lp).(l_p) and .(hp).(h_p) where ll, hh are blocks of binary digits, then lp=lp11l_p = l_{p-1}1 and hp=hp10h_p = h_{p-1}0, and l2p+1=lplp1l_{2p+1} = l_pl_p1 and h2p+1=lphp0h_{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.

  1. input cp1c_{p-1} is a nucleus of period p1p-1 with precision k(p1)k(p-1).

  2. trace external ray lpl_p to depth 3p3p using cp1c_{p-1} as reference.

  3. use Newton’s method to find cpc_p from the endpoint of the ray with precision k(p)k(p).

  4. trace external ray l2p+1l_{2p+1} to depth 3(2p+1)3(2p+1) using cpc_p as reference.

  5. use Newton’s method to find c2p+1c_{2p+1} from the endpoint of the ray with precison k(2p+1)k(2p+1).

  6. set p1:=2p+1p-1 := 2p+1 and repeat from the start.

The function k(p)k(p) to give the required precision from the period may depend on the location, though I suspect it to be O(p)O(p) in general.

The total asymptotic cost of this algorithm is O(pM(k(p)))O(p M(k(p))) which is roughly O(p2.58)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"