# 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…