如何得到无穷小的数字（对于分形）

Question

我正在使用OpenGL用C编程Mandelbrotset，但我遇到了一个问题：我发送到着色器并在着色器中计算的浮点值只能适应一定数量的小数位数。所以如果我放大太远，它就会被像素化。

我曾想过创建自定义数组函数，但我真的搞不懂。除了使用数组，还有其他方法吗？如果不是的话，我怎么能把数组当作一个数字来计算呢？（例如，arr[1,2]x arr[0,2]应给出与仅计算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;
}

********************************************************************************************

Answer 1

您始终可以使用double，但由于这是着色器，将在GPU上执行，因此会带来性能损失。您可以使用的一个技巧是不要使用如此低的值，因为精度会成为一个问题。取而代之的是，当你放大时，你可以在某个点上缩放你的结果，基本上保持数值在稳定的范围内，精度没有问题。IRC，Kerbal空间项目开发人员实际上有一篇关于这项技术的博客文章，所以你可以查看一下。

找不到来自KSP的链接，所以这里有一个类似于Youtube的链接，作者是塞巴斯蒂安·拉格。相关部分约为10分钟标记。

Answer 2

首先，如果您使用双代码而不是浮点代码，您将获得两倍的精度。但是这是GPU代码，所以双代码可以慢一点。（在CPU上，双和浮点的速度差不多）

你要找的东西叫做任意精度算术。要慢得多。您可以找到这样做的库，例如GMP。但是这些是给CPU的。显然，GPU也有一些，我不太熟悉。记住，它很慢！

可以为GPU编写自己的任意精度算法。您不必为每个数组项存储1位数字；您应该存储尽可能多的数据（通常是32位）。

你怎么加？您只需将每个元素加在一起，如果溢出，则带入下一个元素。与减法相同，如果下一个元素下溢，则从该元素借用。

你如何繁殖？就像在小学里，你用每个元素乘以另一个元素，然后在底部把它们加起来。

你如何划分？幸运的是，在Mandelbrot集计算中，没有除法。所以不用麻烦了。

特别是对于Mandelbrot集，有一种更快的近似算法，称为“微扰理论算法”。

微扰理论算法是一种近似计算一组邻近点的Mandelbrot集的方法。首先，选择一个称为“参考点”的点，然后使用通常的算法，以缓慢的任意精度进行计算。参考点必须是黑色像素，或者至少它必须比屏幕上的任何其他像素进行更多的迭代。

公式就在这里——一旦你有了一个参考点，那么对于每个像素，你只需要计算该像素和参考点之间的差值——这是一个很小的数字，所以你不会失去精度！而不是Z

Answer 3

正如其他人建议的那样，使用double而不是float，这将为您提供更高的可能缩放。除此之外，使用分数转义，将允许更高的细节和更少的迭代，因此更大的速度和更好的细节见我的：

• 曼德布罗特

里面有我的Mandelbrot的浮点代码，还有32位和64位浮点的Win32演示链接。然而，双版本着色器不适合回答，所以它们在这里（对于分形，第二遍重新着色着色器并不重要，但你可以从演示中提取它们）：

// 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);
        }
    }

以及：

// 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);
    }

这可以上升到zoom=1e 14，像素开始出现：

在GPU上使用任意精度浮点运算将非常缓慢且有问题（正如其他人已经建议的那样）。然而，有更简单的解决方法来提高浮点或双精度。

例如，您可以将您的值保存为多个双倍的总和。。。

而不是

double x;

您可以使用：

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

其中x=x0-x1-x2。。。xn其中x0保留较小的值，x1较大<代码>xn最大。您只需执行基本的，-，*操作即可

>

x=y

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

x-=y

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

x*=y

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

在每次操作后，您将标准化为每个变量的范围：

   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; }

您需要选择范围，这样您就不会在不使用的数字上浪费变量...

这种精度仍然有限，但精度损失要小得多...

然而，仍然有一个限制（不容易跨越）。现在，您正在从碎片坐标、zoom和panx，y计算位置，这将被限制为浮动，因为我们仍然没有64位双内插器。如果你想打破这个障碍，你需要计算CPU端或顶点上的缩放位置，计算方式与其他计算相同（多个变量之和，但这次是浮动的），并将结果作为变量传递到片段中