I know that quadratic Julia Set matings are only supposed to work for the connected sets, but I got curious and tried to make one using disconnected ones. To my surprise (using the algorithm I learned from this question), the results were really spectacular! My animations are shown here, here, and here through different perspectives (all showing the same mating between the $.285+.01i$ and $-.835-.2321i$ Julia Sets).
Unfortunately, I recently realized why these animations had to be ended early: the orbits for both Julia Sets are unbounded (which is probably why matings only work for certain sets). As such, in the first part of the algorithm when I simply iterate the two coordinates (the $p$ and $q$ sets in the "Initialization" part), they diverge to infinity after about 28 iterations (due to my machine's floating-point precision), which messes up the rest of the calculations. However, the animation clearly shows that the mating could have been continued, and so I'd like to figure out a way to accomplish this.
I realize this could be seen as a programming question as an easy solution could simply be to use arbitrary-precision variables to solve this issue. However, I'd like to know if there are any alternative algorithms to the Thurston Algorithm that avoid these calculations, or perhaps some other magic solution that can accomplish what I seek without quickly diverging to infinity - or at least delay it to a much later iteration.
My project is here to use as testing grounds, if so needed.
Note: I tried implementing robust complex division, but this only saved me a single mating iteration. The orbits diverge much too quickly for this to have any lasting effect.
Edit
I figured I should post what worked for me anyways, even though it uses a programming solution. I made my own BigDouble class, where the mathematical operations performed between numbers are represented in scientific notation. This effectively allows me to bypass the double-precision limit of $1.798*10^{308}$ and go up to $9.99*10^{2,147,483,647}$, all without using slow bignum/infinite-precision arithmetic. Using this, I could now iterate coordinates outside the Mandelbrot Set without running into floating-point limitations (it's not 100% perfect, but it still gets the job done).