# Parallel Mandelbrot Roots

# 1 Serial Algorithm

Newton’s method in one complex variable. Iterate $z$ and $\frac{dz}{dc}$, and update $$c \gets c - \frac{f_c^p(0) - t}{\frac{df_c^p(0)}{dc}}.$$

$O(1)$ space, $O(n)$ work, $O(n)$ time, $O(1)$ processors.

# 2 Matrix Algorithm

Consider the orbit of $z$ from iteration $1$ to iteration $p-1$ as a vector. Then let $f_k = z_k^2 + z_1 - z_{k+1}$ for $k + 1 < p$ and $f_{p-1} = z_{p-1}^2 + z_1 - t$ and solve for the vector $\mathbf{f}(\mathbf{z}) = \mathbf{0}$ using Newton’s method in $p-1$ variables. This involves solving the matrix equation $$\mathbf{J_f} \mathbf{d} = \mathbf{f}$$ and updating $$\mathbf{z} \gets \mathbf{z} - \mathbf{d}.$$

The matrix $\mathbf{J_f}$ has a regular form: it is upper bidiagonal with an added initial column of $1$ values:

$$\mathbf{J_f} = \begin{pmatrix} 1+2z_1 & -1 \\ 1 & 2z_2 & -1 \\ \vdots & & \ddots & \ddots \\ 1 & & & 2z_{n-2} & -1 \\ 1 & & & & 2z_{n-1} \end{pmatrix}$$

We can use a trick based on Sherman-Morrison formula to turn a complicated solve into two simpler solves. As $$\mathbf{J}_f = \mathbf{A} + \begin{pmatrix} 1 & 1 & \cdots & 1 & 1 \end{pmatrix} \otimes \begin{pmatrix} 1 & 0 & \cdots & 0 & 0 \end{pmatrix}$$ so $\mathbf{J_f} \mathbf{d} = \mathbf{f}(\mathbf{z})$ can be solved by solving $$\begin{aligned} \mathbf{A} \mathbf{x} &= \mathbf{f}(\mathbf{z}) \\ \mathbf{A} \mathbf{y} &= \begin{pmatrix} 1 & 1 & \cdots & 1 & 1 \end{pmatrix}^T \\ \mathbf{d} &= \mathbf{x} - \frac{\begin{pmatrix} 1 & 0 & \cdots & 0 & 0 \end{pmatrix}^T \cdot \mathbf{x}}{1 + \begin{pmatrix} 1 & 0 & \cdots & 0 & 0 \end{pmatrix}^T \cdot \mathbf{y}} \mathbf{y} \\ &= \mathbf{x} - \frac{x_1}{1 + y_1} \mathbf{y} \end{aligned}$$

The matrix $\mathbf{A}$ is upper bidiagonal, so can be solved by back-substitution in in $O(n)$ space, $O(n)$ work, $O(n)$ time on $O(1)$ processors.

# 3 Parallel Algorithm

The matrix algorithm can be parallelized to run in in $O(n)$ space, $O(n)$ work, $O(\log(n))$ time on $O(n)$ processors using work-efficient parallel prefix scan.

However, the work done (while still $O(n)$) is a significant constant factor more than the serial algorithm, and the $O(n)$ space requirements make it infeasible for large $n$.

# 4 Benchmark

Task: given the Millionaires #001 deep zoom Mandelbrot set location expressed as a sequence of offsets and periods, find each root along the way to the final view:

1 0.0000000000000000e+00 0.0000000000000000e+00
2 -1.0000000000000000e+00 0.0000000000000000e+00
3 -7.5487766624669272e-01 0.0000000000000000e+00
4 -1.8592214028279200e-01 0.0000000000000000e+00
8 -9.8228378871340798e-04 0.0000000000000000e+00
69 2.1724309435420009e-04 2.3489118850327509e-04
142 1.3225292985888009e-11 7.1358633969614633e-12
272 2.8529310483870992e-16 1.0267761000262062e-15
560 -1.6383387496504541e-22 2.0300116425226282e-22
1104 -2.1353430107983675e-33 7.6284717085111253e-34
2228 1.5911948904690584e-49 4.0112501481921464e-49
4436 2.3215027144370148e-74 1.3325698720132981e-73
8896 -5.9289281961504007e-110 9.1490695754906453e-110
17768 -1.6130749406127237e-164 3.3282600257703784e-165
35564 7.5630380778794724e-247 2.4419891481973179e-246
71100 -2.2206801493397202e-370 1.8281359915497428e-369
142232 -4.3057238712082329e-554 3.0489078189532703e-554
284432 4.8325042482703457e-831 2.4686705119967285e-831
568900 -1.1303349232678663e-1246 -9.4216823889086909e-1247
1137764 -6.2022423569065681e-1870 5.0076670707505569e-1870
5.0e-2805

The final line is the size of the view; the precision required for each root is based on the magnitude of the next line’s offset.

The serial algorithm is trivial to implement, and took 5 minutes at 99% CPU with 4MB peak RAM usage.

The parallel algorithm is hard to implement, and took 6min40s at 1548% CPU with 17GB peak RAM usage.

Conclusion: the serial algorithm is better in all aspects, at least on a CPU with modest parallelism (16 threads). On a GPU the story might be different. Work in progress…