Optimizing Complex Mobius Transformations on a Fragment Shader

Question

I'm developing my own graphics engine to render all sorts of fractals (like my video here for example), and I'm currently working on optimizing my code for Julia Set matings (see this question and my project for more details). In the fragment shader, I use this function:

vec3 JuliaMatingLoop(dvec2 z)
{
    ...

    for (int k = some_n; k >= 0; --k)
    {
        // z = z^2
        z = dcproj(c_2(z));
        
        // Mobius Transformation: (ma[k] * z + mb[k]) / (mc[k] * z + md[k])
        z = dcproj(dc_div(cd_mult(ma[k], z) + mb[k], dc_mult(mc[k], z) + md[k]));
    }
    
    ...
}

And after reading this, I realized that I'm doing Mobius transformations in this code, and (mathematically speaking) I can use matrices to accomplish the same operation. However, the a, b, c, and d constants are all complex numbers (represented as ma[k], mb[k], mc[k], and md[k] in my code), whereas the elements in GLSL matrices contain only real numbers (rather than vec2). And so to my question: is there a way to optimize these Mobius transformations using matrices in GLSL? Or any other way of optimizing this part of my code?

Helper functions (I need to use doubles for this part, so I can't optimize by switching to using floats):

// Squaring
dvec2 c_2(dvec2 c)
{
    return dvec2(c.x*c.x - c.y*c.y, 2*c.x*c.y);
}
// Multiplying
dvec2 dc_mult(dvec2 a, dvec2 b)
{
    return dvec2(a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x);
}
// Dividing
dvec2 dc_div(dvec2 a, dvec2 b)
{
    double x = b.x * b.x + b.y * b.y;
    return vec2((a.x * b.x + a.y * b.y) / x, (b.x * a.y - a.x * b.y) / x);
    
}
// Riemann Projecting
dvec2 dcproj(dvec2 c)
{
    if (!isinf(c.x) && !isinf(c.y) && !isnan(c.x) && !isnan(c.y))
        return c;
        
    return dvec2(infinity, 0);
}

Answer 1

I'm not sure if this will help, but yes you can do complex arithmetic by matrices.

If you regard a complex number z as a real two-vector with components Re(z), Im(z) Then

A*z + B ~ (Re(A)  -Im(A) ) * (Re(z)) + (Re(B))
          (Im(A)   Re(A) )   (Im(z))   (Im(B))

Of course you actually want

(A*z + B) / (C*z + D)

If you compute

A*z+b as (x)
         (y)

C*z+d as (x')
         (y')

Then the answer you seek is

inv( x' -y') * ( x)
   ( y' x' )   ( y)
i.e
    (1/(x'*x'+y'*y')) * (x'  y') * (x)
                        (-y' x')   (y)

One thing to note, though, is that in these formulae, as in your code, division is not implemented as robustly as it could be. The trouble lies in evaluating b.x * b.x + b.y * b.y. This could overflow to infinity, or underflow to 0, even though the result of division could be quite reasonable. A commonly used way round this is Smith's method eg here and if you search for 'robust complex division' you'll find more recent work. Often this sort of thing matters little, but if you are iterating off to infinity it could make a difference.