Algorithm for computing external angles for the Mandelbrot set

Question

Let $M$ be the Mandelbrot set: there exists a unique series $$ \psi(z) := z + \sum_{m=0}^{+\infty} b_m z^{-m} = z - \frac{1}{2} + \frac{1}{8} z^{-1} - \frac{1}{4} z^{-2} + \cdots $$ which defines a conformal bijection between the complement $\{z\in\mathbb{C} : |z|>1\}$ of the closed unit disk in $\mathbb{C}$ and the complement of the Mandelbrot set. This function $\psi$ is sometimes known as the “Jungreis function”, see the answers to this question for more. The argument $\arg(\psi^{-1}(w))$ for a point $w\not\in M$ is called the external angle of $w$.

There exists an easy way to “almost” compute $\arg(\psi^{-1}(w))$, or $\psi^{-1}(w)$ itself for that matter: indeed, if we let $p_0(w) := w$ and $p_{i+1}(w) := p_i(w)^2 + w$, then $$ p_n(\psi(z)) = z^{2^n} + o(1) $$ as $z\to\infty$, so $\psi^{-1}(w)$ can be “almost” computed as the limit of the $(2^n)$-th root of $p_n(w)$ as $n\to+\infty$. The reason for the “almost” is that, while this indeed allows for computation of the modulus $|\psi^{-1}(w)|$, it leads to an indetermination between $2^n$ values on the argument. The formula $$ \psi^{-1}(w) = w \mskip3mu \prod_{n=1}^{+\infty} \left(1 + \frac{w}{p_{n-1}(w)^2}\right)^{1/2^n} $$ is no better (it also requires computing $(2^n)$-th roots and one cannot simply take the principal determination; my understanding is that one needs to find a determination of $\Big(1 + \frac{w}{p_{n-1}(w)^2}\Big)^{1/2^n}$ that is continuous outside of $M$, which seems computationally intractable).

So, is there a way to lift this square root indeterminacy and compute external arguments for arbitrary $w\not\in M$? Is there an algorithm that does this in a reasonably efficient way (which excludes, e.g., trying to trace external rays outwards towards infinity)?

I was unable to find anything relevant in the literature. There is a 1986 paper by Douady titled “Algorithms for computing angles in the Mandelbrot set” which seems promising, but it seems to concerns the computation for points of $M$, not points outside $M$. This web page about the Mandel program actually discusses the issue (in the section called “Computation of the external argument”), but the description is vague (e.g., where it speaks of a “modified” function $\arg(z/(z-c))$), and the conclusion that “the discontinuities are moved closer to the Mandelbrot set” is not too promising.

Answer 1

I have implemented some algorithms based on Kawahira's paper, which as presented goes $\theta \to c \not\in M$, but can be adapted to go $c \to \theta$. $\theta$ is conveniently expressed in turns as a binary expansion. When tracing inwards, one peels off the most-significant bit (aka angle doubling) each time the ray crosses a dwell band (integer part of normalized iteration count increases by 1). The trick when tracing outwards is to prepend bits when crossing dwell bands, depending if the outer cell was entered from its left or right inner cell. A picture may make it clearer:

The exterior grid is generated from the fractional part of the smoothed iteration count, and the argument of the final (first to escape) $z$ iterate. One can see that approaching the $\frac{1}{2}$ bond point, the cells alternate left/right corresponding to the binary expansion $.(01)$ or $.(10)$. The argument of the first iterate to escape (typically floating point), within the cell of the starting point, can be used to get a few more least significant bits for the accumulated angle, but the prefix is found by accumulating bits one by one when tracing the ray outwards.

Source code:

• m_d_exray_in.c $\theta \to c$, machine double precision

• m_r_exray_in.c $\theta \to c$ in arbitrary (dynamically changed as necessary) precision

• m_d_exray_out.c $c \to \theta$, machine double precision

• m_r_exray_out.c $c \to \theta$, arbitrary (but fixed, dynamic is possible but still TODO) precision

• m-exray-in.c command line driver program showing usage of the library functions

• m-exray-out.c command line driver program showing usage of the library functions

git clone https://code.mathr.co.uk/mandelbrot-numerics.git

However, tracing external rays to/from dwell $n$ takes $O(n^2)$ time (even ignoring that higher $n$ needs higher precision which costs more), which may make it too slow in practice. I am also looking for faster methods since some years, but I haven't found any yet. See my related question: fast algorithms for external angle computations

Answer 2

Googling for the paper of Douady, I found this paper by Tomoki Kawahira: he explicitly computes the quantity $\psi^{-1}(w)$ for each $w\in\Bbb C\setminus M$ by using Newton's method (he considers $\Phi=\psi^{-1}$), and this allows also an error estimate and a quite fast convergence rate: could it be of some help?