how do I get infinitely small numbers (for fractals)

Question

I'm programming the Mandelbrotset with C++ using OpenGL, but I've run into a problem: the floats I'm sending to and calculating in my shader can only fit a certain amount of decimal places. So if I zoom in too far it just gets pixelated.

I thought about creating custom array functions but I can't really figure it out. Is there any other way than using arrays? and if not how can I calculate stuff using arrays as if they are a single number? (e.g. arr[1,2] x arr[0,2] should give the same output as just calculating 1.2 x 0.2)

in vec4 gl_FragCoord;
 
out vec4 frag_color;
uniform float zoom;
uniform float x;
uniform float y;
#define maximum_iterations 1000

int mandelbrot_iterations()
{
    float real = ((gl_FragCoord.x / 1440.0f - 0.5f) * zoom + x) * 4.0f;
    float imag = ((gl_FragCoord.y / 1440.0f - 0.5f) * zoom + y) * 4.0f;
 
    int iterations = 0;
    float const_real = real;
    float const_imag = imag;
 
    
    while (iterations < maximum_iterations)
    {
        float tmp_real = real;
        real = (real * real - imag * imag) + const_real;
        imag = (2.0f * tmp_real * imag) + const_imag;
         
        float dist = real * real + imag * imag;
 
        if (dist > 4.0f)
            break;
 
        ++iterations;
    }
    return iterations;
}

^the Mandelbrot function in my fragmentshader

Marcus

You could start by using double instead of float, that gives you double precision. If you need more than that, you can look up 'arbitrary precision arithmetic' and find a library like GMP. But be warned: arbitrary precision arithmetic is sloooooow. You can also look up the 'perturbation theory approach' to calculating the Mandelbrot set. – Coomb

@Coomb doubles aren't supported by OpenGL, but the arbitrary precision arithmetic looks promising and I found a GLSL library for it but I don't really understand how it works and speed is important as I'm trying to show it realtime (so I can show it on screen in my school presentation). – Marcus

@SpaceKek double is supported by OpenGL (you need to use core profile or enable extention for older GLSL) see my answer using double precision however you're partialy right doubles are not implemented for interpolators so varying variables passed between shaders can not be double... – Enswathe

Answer 1

As others suggested use double instead of float that will give you much higher possible zoom. On top of that use fractional escape that will allow much higher detail with much less iterations hence bigger speed with better details see mine:

• GLSL Mandelbrot

In there is float code for mine Mandelbrot And links to Win32 demos for both 32 and 64 bit float. However double version shaders did not fit to answer so here they are (for the fractal, the second pass recoloring shaders are not important but you can extract them form the demo):

// Fragment
#version 450 core
uniform dvec2 p0=vec2(0.0,0.0);     // mouse position <-1,+1>
uniform double zoom=1.000;          // zoom [-]
uniform int  n=100;                 // iterations [-]
uniform int  sh=7;                  // fixed point accuracy [bits]
uniform int  multipass=0;           // multi pass?
uniform int  inverted=0;            // inverted/reciprocal position?
in smooth vec2 p32;
out vec4 col;

const int n0=1;                     // forced iterations after escape to improve precision

vec3 spectral_color(float l)        // RGB <0,1> <- lambda l <400,700> [nm]
    {
    float t;  vec3 c=vec3(0.0,0.0,0.0);
         if ((l>=400.0)&&(l<410.0)) { t=(l-400.0)/(410.0-400.0); c.r=    +(0.33*t)-(0.20*t*t); }
    else if ((l>=410.0)&&(l<475.0)) { t=(l-410.0)/(475.0-410.0); c.r=0.14         -(0.13*t*t); }
    else if ((l>=545.0)&&(l<595.0)) { t=(l-545.0)/(595.0-545.0); c.r=    +(1.98*t)-(     t*t); }
    else if ((l>=595.0)&&(l<650.0)) { t=(l-595.0)/(650.0-595.0); c.r=0.98+(0.06*t)-(0.40*t*t); }
    else if ((l>=650.0)&&(l<700.0)) { t=(l-650.0)/(700.0-650.0); c.r=0.65-(0.84*t)+(0.20*t*t); }
         if ((l>=415.0)&&(l<475.0)) { t=(l-415.0)/(475.0-415.0); c.g=             +(0.80*t*t); }
    else if ((l>=475.0)&&(l<590.0)) { t=(l-475.0)/(590.0-475.0); c.g=0.8 +(0.76*t)-(0.80*t*t); }
    else if ((l>=585.0)&&(l<639.0)) { t=(l-585.0)/(639.0-585.0); c.g=0.84-(0.84*t)           ; }
         if ((l>=400.0)&&(l<475.0)) { t=(l-400.0)/(475.0-400.0); c.b=    +(2.20*t)-(1.50*t*t); }
    else if ((l>=475.0)&&(l<560.0)) { t=(l-475.0)/(560.0-475.0); c.b=0.7 -(     t)+(0.30*t*t); }
    return c;
    }

void main()
    {
    int i,j,N;
    dvec2 pp,p;
    double x,y,q,xx,yy,mu,cx,cy;
    p=dvec2(p32);

    pp=(p/zoom)-p0;         // y (-1.0, 1.0)
    pp.x-=0.5;              // x (-1.5, 0.5)
    if (inverted!=0)
        {
        cx=pp.x/((pp.x*pp.x)+(pp.y*pp.y));  // inverted
        cy=pp.y/((pp.x*pp.x)+(pp.y*pp.y));
        }
    else{
        cx=pp.x;                // normal
        cy=pp.y;
        }
    for (x=0.0,y=0.0,xx=0.0,yy=0.0,i=0;(i<n-n0)&&(xx+yy<4.0);i++)
        {
        q=xx-yy+cx;
        y=(2.0*x*y)+cy;
        x=q;
        xx=x*x;
        yy=y*y;     
        }
    for (j=0;j<n0;j++,i++)  // 2 more iterations to diminish fraction escape error
        {
        q=xx-yy+cx;
        y=(2.0*x*y)+cy;
        x=q;
        xx=x*x;
        yy=y*y;
        }
    mu=double(i)-double(log2(log(float(sqrt(xx+yy)))));
    mu*=double(1<<sh); i=int(mu);
    N=n<<sh;
    if (i>N) i=N;
    if (i<0) i=0;

    if (multipass!=0)
        {
        // i
        float r,g,b;
        r= i     &255; r/=255.0;
        g=(i>> 8)&255; g/=255.0;
        b=(i>>16)&255; b/=255.0;
        col=vec4(r,g,b,255);
        }
    else{
        // RGB
        float q=float(i)/float(N);
        q=pow(q,0.2);
        col=vec4(spectral_color(400.0+(300.0*q)),1.0);
        }
    }

And:

// Vertex
#version 450 core
layout(location=0) in vec2 pos;         // glVertex2f <-1,+1>
out smooth vec2 p32;                    // texture end point <0,1>
void main()
    {
    p32=pos;
    gl_Position=vec4(pos,0.0,1.0);
    }

This can go up to zoom = 1e+14 where pixelations starts to occur:

Using arbitrary precision float on GPU will be very slow and problematic (as others already suggested). However there are simpler workarounds to enhance precision of float or double.

For example you can hold your value as sum of more doubles...

instead of

double x;

you can use:

double x0,x1,x2,....,xn;

where x = x0+x1+x2+...+xn where x0 hold small values, x1 bigger , ... xn biggest. You need just basic +,-,* operations so

1. x += y

   x0+=y0; x1+=y1; ... xn+=yn;
   

2. x -= y

   x0-=y0; x1-=y1; ... xn-=yn;
   

3. x *= y

   x0*=(y0+y1+...yn);
   x1*=(y0+y1+...yn);
   ...
   xn*=(y0+y1+...yn);
   

And after every operation you normalize to ranges of each variable:

   if (fabs(x0)>1e-10){ x1+=x0; x0=0; }
   if (fabs(x1)>1e-9) { x2+=x1; x1=0; }
   ...
   if (fabs(x(n-1))>1e+9){ xn+=x(n-1); x(n-1)=0; }

You need to chose the ranges so you not waste variables on numbers that will not be used ...

With this precision is still limited but the accuracy loss is much much smaller...

However there is still one limit (which is not crossable easily). Right now you are computing position from fragment coordinate, zoom and pan x,y which will be limited to float as we still do not have 64bit double interpolators. If you want break this barrier you need to compute the zoomed position either on CPU side or on Vertex in the same manner as rest of the computations (sum of more variables but float this time) and pass the result as varying into fragment

Enswathe

Hi, I tried to reach you by email couldn't, would you help me please with that question :) #65025929 – Polonium

Answer 2

Well first, if you use double instead of float, you get twice as much precision. But this is GPU code, so double can be slower. (On a CPU, double and float are about the same speed)

The thing you're looking for is called arbitrary precision arithmetic. It's a lot slower. You can find libraries which do this, such as GMP. But those are for the CPU. Apparently, there are also some for the GPU - I'm not familiar with them. Remember, it's slow!

You could write your own arbitrary precision arithmetic for the GPU, yes. You don't have to store 1 digit per array item; you should store as many as possible (usually 32 bits).

How do you add? You just add together each element, and carry into the next element if it overflows. Same for subtraction, you borrow from the next element if it underflows.

How do you multiply? Like in grade school - you multiply each element by each other element, then you add them together at the bottom.

How do you divide? Luckily, in the Mandelbrot set calculation, there is no division. So don't bother.

---

For the Mandelbrot set specifically, there is a faster approximate algorithm which is just known as "the perturbation theory algorithm".

The perturbation theory algorithm is an approach to approximately calculate the Mandelbrot set for a bunch of nearby points. First, you pick one point called the reference point, and you calculate it using the usual algorithm, with slow arbitrary precision maths. The reference point must be a black pixel, or at least it must go for more iterations than any other pixel on the screen.

The formula is here - once you have a reference point, then for each pixel, you only need to calculate the difference between that pixel and the reference point - which is a small number, so you don't lose precision! Instead of Z <- Z^2 + C you calculate Zref+Zdiff <- (Zref+Zdiff)^2 + (Cref+Cdiff) - if you simplify this equation mathematically, you will find that a lot of the Zref's and Cref's cancel out, and you can calculate Zdiff without losing precision. (Don't just calculate it the normal way and then subtract Zref, because you won't get any extra precision that way)

Coomb

Answer 3

You can always go with double, but since this is shader and will be executed on GPU it will come with a performance penalty. One trick you can use is to not go to such low values, where precision becomes a problem. Instead as you zoom in you can at some point scale your results back up, basically keeping the values in steady range where precision is not a problem. IRC, Kerbal Space Program developers actually have a blog post about this technique, so you can check it out.

Can't find the one from KSP, so here is a link to something similar on Youtube, by Sebastian Lague. The relevant part is around 10 minute mark.

Iy

This doesn't fix the problem. The Mandelbrot set calculation requires the use of large numbers and small numbers together. You still lose precision by doing math with very large numbers and "normal-sized" numbers together. – Coomb