Search code examples
c++openglglslfractalsmandelbrot

how do I get infinitely small numbers (for fractals)

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

Solution

  • 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:

    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:

    pixelation

    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