我如何获得无限小的数字（对于分形）
how do I get infinitely small numbers (for fractals)
我正在使用 C++ 使用 OpenGL 对 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;
}
^我的片段着色器中的 Mandelbrot 函数
您始终可以使用 double，但由于这是 shader 并且将在 GPU 上执行，因此会带来性能损失。您可以使用的一个技巧是不要使用如此低的值，否则精度会成为问题。相反，当您放大时，您可以在某些时候将结果放大，基本上将值保持在精度不成问题的稳定范围内。 IRC，Kerbal Space 程序开发人员实际上有一篇关于此技术的博客 post，因此您可以查看。
无法从 KSP 找到，所以这里是 Youtube 上的 link 到 something similar，作者是 Sebastian Lague。相关部分大约是 10 分钟标记。
首先，如果您使用 double 而不是 float，您将获得两倍的精度。但这是 GPU 代码，因此 double 可能会更慢。 （在CPU上，double和float速度差不多）
你要找的东西叫做任意精度算术。它慢了很多。您可以找到执行此操作的库，例如 GMP。但那些是给 CPU 的。显然，还有一些用于 GPU - 我对它们不熟悉。记住，它很慢！
您可以为 GPU 编写您自己的任意精度算法，是的。您不必为每个数组项存储 1 个数字；你应该存储尽可能多的（通常是 32 位）。
如何添加？您只需将每个元素加在一起，如果它溢出，则将其带入下一个元素。减法也一样，如果下溢，则从下一个元素借用。
你如何乘法？就像在小学一样 - 你将每个元素乘以其他元素，然后在底部将它们相加。
你是怎么分的？幸运的是，在 Mandelbrot 集合计算中，没有除法。所以不要打扰。
具体针对Mandelbrot集，有一种更快的近似算法，即所谓的“微扰理论算法”。
微扰理论算法是一种近似计算一组附近点的Mandelbrot集的方法。首先，你选择一个称为参考点的点，然后使用通常的算法计算它，使用缓慢的任意精度数学。参考点必须是黑色像素，或者至少它必须比屏幕上的任何其他像素进行更多次迭代。
公式是here - 一旦你有了一个参考点，那么对于每个像素，你只需要计算该像素与参考点之间的差异 - 这是一个很小的数字，所以你不需要'失去精度！您计算的不是 Z <- Z^2 + C，而是 Zref+Zdiff <- (Zref+Zdiff)^2 + (Cref+Cdiff) - 如果您在数学上简化此等式，您会发现很多 Zref 和 Cref 抵消了，并且您可以在不损失精度的情况下计算 Zdiff。 （不要只是用正常的方式计算它然后减去 Zref，因为那样你不会得到任何额外的精度）
正如其他人建议的那样，使用 double 而不是 float，这会给您带来更高的缩放比例。最重要的是使用分数转义，它将允许更多的细节和更少的迭代因此更快的速度和更好的细节见我的：
其中有用于我的 Mandelbrot 的浮点代码以及指向 32 位和 64 位浮点的 Win32 演示的链接。但是 double 版本着色器不适合回答，所以在这里（对于分形，第二遍重新着色着色器并不重要，但您可以从演示中提取它们）：
// 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 上使用任意精度的浮点数会非常缓慢且有问题（正如其他人已经建议的那样）。然而，有更简单的解决方法来提高 float 或 double 的精度。
例如，您可以将您的价值保存为更多 doubles...
的总和
而不是
double x;
您可以使用：
double x0,x1,x2,....,xn;
where x = x0+x1+x2+...+xn where x0 hold small values, x1 bigger , ... xn maximum.您只需要基本的 +,-,* 操作，所以
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 和平移 x,y 计算位置，这将被限制为浮动，因为我们仍然没有 64 位双插值器。如果你想打破这个障碍，你需要以与其余计算相同的方式计算 CPU 侧或 Vertex 上的缩放位置（更多变量的总和，但这次是浮动的）并将结果作为变化传递给片段
我正在使用 C++ 使用 OpenGL 对 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;
}
^我的片段着色器中的 Mandelbrot 函数
您始终可以使用 double，但由于这是 shader 并且将在 GPU 上执行，因此会带来性能损失。您可以使用的一个技巧是不要使用如此低的值，否则精度会成为问题。相反，当您放大时，您可以在某些时候将结果放大，基本上将值保持在精度不成问题的稳定范围内。 IRC，Kerbal Space 程序开发人员实际上有一篇关于此技术的博客 post，因此您可以查看。
无法从 KSP 找到，所以这里是 Youtube 上的 link 到 something similar，作者是 Sebastian Lague。相关部分大约是 10 分钟标记。
首先，如果您使用 double 而不是 float，您将获得两倍的精度。但这是 GPU 代码，因此 double 可能会更慢。 （在CPU上，double和float速度差不多）
你要找的东西叫做任意精度算术。它慢了很多。您可以找到执行此操作的库，例如 GMP。但那些是给 CPU 的。显然，还有一些用于 GPU - 我对它们不熟悉。记住，它很慢！
您可以为 GPU 编写您自己的任意精度算法，是的。您不必为每个数组项存储 1 个数字；你应该存储尽可能多的（通常是 32 位）。
如何添加？您只需将每个元素加在一起，如果它溢出，则将其带入下一个元素。减法也一样，如果下溢，则从下一个元素借用。
你如何乘法？就像在小学一样 - 你将每个元素乘以其他元素，然后在底部将它们相加。
你是怎么分的？幸运的是，在 Mandelbrot 集合计算中，没有除法。所以不要打扰。
具体针对Mandelbrot集，有一种更快的近似算法，即所谓的“微扰理论算法”。
微扰理论算法是一种近似计算一组附近点的Mandelbrot集的方法。首先，你选择一个称为参考点的点，然后使用通常的算法计算它，使用缓慢的任意精度数学。参考点必须是黑色像素，或者至少它必须比屏幕上的任何其他像素进行更多次迭代。
公式是here - 一旦你有了一个参考点，那么对于每个像素，你只需要计算该像素与参考点之间的差异 - 这是一个很小的数字，所以你不需要'失去精度！您计算的不是 Z <- Z^2 + C，而是 Zref+Zdiff <- (Zref+Zdiff)^2 + (Cref+Cdiff) - 如果您在数学上简化此等式，您会发现很多 Zref 和 Cref 抵消了，并且您可以在不损失精度的情况下计算 Zdiff。 （不要只是用正常的方式计算它然后减去 Zref，因为那样你不会得到任何额外的精度）
正如其他人建议的那样，使用 double 而不是 float，这会给您带来更高的缩放比例。最重要的是使用分数转义，它将允许更多的细节和更少的迭代因此更快的速度和更好的细节见我的：
其中有用于我的 Mandelbrot 的浮点代码以及指向 32 位和 64 位浮点的 Win32 演示的链接。但是 double 版本着色器不适合回答，所以在这里（对于分形，第二遍重新着色着色器并不重要，但您可以从演示中提取它们）：
// 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 上使用任意精度的浮点数会非常缓慢且有问题（正如其他人已经建议的那样）。然而，有更简单的解决方法来提高 float 或 double 的精度。
例如，您可以将您的价值保存为更多 doubles...
而不是
double x;
您可以使用：
double x0,x1,x2,....,xn;
where x = x0+x1+x2+...+xn where x0 hold small values, x1 bigger , ... xn maximum.您只需要基本的 +,-,* 操作，所以
x += yx0+=y0; x1+=y1; ... xn+=yn;x -= yx0-=y0; x1-=y1; ... xn-=yn;x *= yx0*=(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 和平移 x,y 计算位置，这将被限制为浮动，因为我们仍然没有 64 位双插值器。如果你想打破这个障碍，你需要以与其余计算相同的方式计算 CPU 侧或 Vertex 上的缩放位置（更多变量的总和，但这次是浮动的）并将结果作为变化传递给片段