渲染 mandelbrot 集的最快算法是什么？
What are the fastest algorithms for rendering the mandelbrot set?
我已经尝试了很多算法来渲染 Mandelbrot 集，包括朴素的逃逸时间算法，以及优化的逃逸时间算法。但是，是否有更快的算法可以像我们在 YouTube 上看到的那样有效地产生真正的深度缩放。此外，我很想获得一些关于如何提高 C/C++ double
精度的想法
优化后的转义算法应该足够快，可以实时绘制 Mandelbrot 集。您可以使用 多线程 以便您的实施速度更快（例如，使用 OpenMP 非常容易）。如果需要，您还可以使用 SIMD 指令 手动矢量化您的代码以使其更快。如果这对您来说仍然不够快（尽管要有效地做到这一点有点复杂）。最后，您应该调整迭代次数，使其相当小。
缩放不应该对性能有任何直接影响。它只是改变了计算的输入window。但是，它确实具有 间接 影响，因为实际迭代次数会发生变化。 不应计算 window 之外的点。
双精度也应该足以正确绘制 Mandelbrot 集。但是如果你真的想要更精确的计算，你可以使用double-double precision，它提供了相当好的精度并且性能也不错。但是，手动实现 double-double precision 有点棘手，而且它仍然比仅使用 double precision 要慢得多。
我最快的解决方案通过遵循等高线边界和填充来避免在相同深度的大面积上进行迭代。有一个惩罚是可以剪掉小芽而不是绕过它们，但总而言之，为快速变焦付出的代价很小。
一个可能的效率是，如果缩放比例加倍，您已经拥有 ¼ 的点。
对于动画，我记录了每一帧的值，每次将比例加倍，并在播放时实时插入中间帧，因此动画每秒加倍一次。 double 类型允许存储超过 50 个关键帧，给出一个持续超过一分钟的动画（进入然后退出）。
实际迭代是由手工汇编程序完成的，因此一个像素完全在 FPU 中迭代。
即使是高端 CPU 与普通 GPU 相比也会慢得多。即使在 GPU 上使用简单的迭代算法，您也可以获得实时渲染。因此，在 GPU 上使用更好的算法可以获得高变焦，但是对于您需要的任何体面的算法：
- 多通道渲染，因为我们不能在 GPU 上自行修改纹理
- 高精度浮点数
floats/doubles 是不够的。
这里有几个相关的问答：
- Interior Distance Estimate algorithm for the Mandelbrot set
- Perturbation theory
这可能会让你开始......
一种加快速度的方法是您可以像我在第一个 link 中那样使用 分数转义 。它提高了图像质量，同时保持较低的最大迭代次数。
第二个 link 将让您近似分形的哪些部分进出以及多远。它不是很准确，但可以用来避免计算“肯定在外面”的部分的迭代。
接下来 link 将向您展示如何获得更高的精度。
最后一个 link 是关于 扰动 这个想法是你只对一些参考点使用高精度数学，并使用它来计算它的低精度数学的邻居点在不失去精度的情况下。从未使用过，但看起来很有前途。
最后，一旦您实现了快速渲染，您可能希望以此为目标：
- How to adjust panning while zooming Mandelbrot set
这里是 GLSL 中用于单值的 3* 64 位双精度的小例子：
// high precision float (very slow)
dvec3 fnor(dvec3 a)
{
dvec3 c=a;
if (abs(c.x)>1e-5){ c.y+=c.x; c.x=0.0; }
if (abs(c.y)>1e+5){ c.z+=c.y; c.y=0.0; }
return c;
}
double fget(dvec3 a){ return a.x+a.y+a.z; }
dvec3 fset(double a){ return fnor(dvec3(a,0.0,0.0)); }
dvec3 fadd(dvec3 a,double b){ return fnor(a+fset(b)); }
dvec3 fsub(dvec3 a,double b){ return fnor(a-fset(b)); }
dvec3 fmul(dvec3 a,double b){ return fnor(a*b); }
dvec3 fadd(dvec3 a,dvec3 b){ return fnor(a+b); }
dvec3 fsub(dvec3 a,dvec3 b){ return fnor(a-b); }
dvec3 fmul(dvec3 a,dvec3 b)
{
dvec3 c;
c =fnor(a*b.x);
c+=fnor(a*b.y);
c+=fnor(a*b.z);
return fnor(c);
}
所以每个高精度值都是 dvec3 ... fnor 中的阈值可以更改为任何范围。您可以将其转换为 vec3 和 float ...
[Edit1]“快速”C++ 示例
好的，我想试试我的新 along with my AVR32 MCU to compute Mandelbrot So I can compare speed with this Arduino + 3D + Pong + Mandelbrot。我使用 AT32UC3A3256 ~66MHz 无 FPU 无 GPU 和 128x64x1bpp 显示器。没有外部存储器只有内部 16+32+32 KByte。 Naive Mandlebrot 的速度很慢（每帧约 2.5 秒）所以我搞砸了这样的事情（利用那个位置和视图的缩放是连续的）：
将分辨率降低 2
为抖动腾出空间，因为我的输出只是黑白
基于zoom
使用变量最大迭代n
更改 n 后会使最后一帧无效以强制执行完全重新计算。我知道这很慢，但在缩放范围之间的转换中只发生 3 次。
最后一帧的缩放计数看起来不太好，因为它不是线性的。
可以使用最后的计数，但为此还需要记住用于迭代的复杂变量，这会占用太多内存。
记住最后一帧以及 x,y 屏幕坐标映射到哪个 Mandelbrot 坐标。
在每一帧计算屏幕坐标和 Mandelbrot 坐标之间的映射。
重新映射最后一帧以调整到新位置和缩放
所以只需查看来自 #3、#4 的数据，如果我们在最后一帧和实际帧中的位置相同（接近像素大小的一半），复制像素。并重新计算其余部分。
如果您的视图平滑（因此位置和缩放在每帧基础上不会有太大变化），这将极大地提高性能。
我知道它的描述有点模糊，所以这里有一个C++代码，您可以在其中推断所有疑问：
//---------------------------------------------------------------------------
//--- Fast Mandelbrot set ver: 1.000 ----------------------------------------
//---------------------------------------------------------------------------
template<int xs,int ys,int sh> void mandelbrot_draw(float mx,float my,float zoom)
{
// xs,ys - screen resolution
// sh - log2(pixel_size) ... dithering pixel size
// mx,my - Mandelbrot position (center of view) <-1.5,+0.5>,<-1.0,+1.0>
// zoom - zoom
// ----------------
// (previous/actual) frame
static U8 p[xs>>sh][ys>>sh]; // intensities (raw Mandelbrot image)
static int n0=0; // max iteraions
static float px[(xs>>sh)+1]={-1000.0}; // pixel x position in Mandlebrot
static float py[(ys>>sh)+1]; // pixel y position in Mandlebrot
// temp variables
U8 shd; // just pattern for dithering
int ix,iy,i,n,jx,jy,kx,ky,sz; // index variables
int nx=xs>>sh,ny=ys>>sh; // real Mandelbrot resolution
float fx,fy,fd; // floating Mandlebrot position and pixel step
float x,y,xx,yy,q; // Mandelbrot iteration stuff (this need to be high precision)
int qx[xs>>sh],qy[ys>>sh]; // maping of pixels between last and actual frame
float px0[xs>>sh],py0[ys>>sh]; // pixel position in Mandlebrot from last frame
// init vars
if (zoom< 10.0) n= 31;
else if (zoom< 100.0) n= 63;
else if (zoom< 1000.0) n=127;
else n=255;
sz=1<<sh;
ix=xs; if (ix>ys) ix=ys; ix/=sz;
fd=2.0/(float(ix-1)*zoom);
mx-=float(xs>>(1+sh))*fd;
my-=float(ys>>(1+sh))*fd;
// init buffers
if ((px[0]<-999.0)||(n0!=n))
{
n0=n;
for (ix=0;ix<nx;ix++) px[ix]=-999.0;
for (iy=0;iy<ny;iy++) py[iy]=-999.0;
for (ix=0;ix<nx;ix++)
for (iy=0;iy<ny;iy++)
p[ix][iy]=0;
}
// store old and compute new float positions of pixels in Mandelbrot to px[],py[],px0[],py0[]
for (fx=mx,ix=0;ix<nx;ix++,fx+=fd){ px0[ix]=px[ix]; px[ix]=fx; qx[ix]=-1; }
for (fy=my,iy=0;iy<ny;iy++,fy+=fd){ py0[iy]=py[iy]; py[iy]=fy; qy[iy]=-1; }
// match old and new x coordinates to qx[]
for (ix=0,jx=0;(ix<nx)&&(jx<nx);)
{
x=px[ix]; y=px0[jx];
xx=(x-y)/fd; if (xx<0.0) xx=-xx;
if (xx<=0.5){ qx[ix]=jx; px[ix]=y; }
if (x<y) ix++; else jx++;
}
// match old and new y coordinates to qy[]
for (ix=0,jx=0;(ix<ny)&&(jx<ny);)
{
x=py[ix]; y=py0[jx];
xx=(x-y)/fd; if (xx<0.0) xx=-xx;
if (xx<=0.5){ qy[ix]=jx; py[ix]=y; }
if (x<y) ix++; else jx++;
}
// remap p[][] by qx[]
for (ix=0,jx=nx-1;ix<nx;ix++,jx--)
{
i=qx[ix]; if ((i>=0)&&(i>=ix)) for (iy=0;iy<ny;iy++) p[ix][iy]=p[i][iy];
i=qx[jx]; if ((i>=0)&&(i<=jx)) for (iy=0;iy<ny;iy++) p[jx][iy]=p[i][iy];
}
// remap p[][] by qy[]
for (iy=0,jy=ny-1;iy<ny;iy++,jy--)
{
i=qy[iy]; if ((i>=0)&&(i>=iy)) for (ix=0;ix<nx;ix++) p[ix][iy]=p[ix][i];
i=qy[jy]; if ((i>=0)&&(i<=jy)) for (ix=0;ix<nx;ix++) p[ix][jy]=p[ix][i];
}
// Mandelbrot
for (iy=0,ky=0,fy=py[iy];iy<ny;iy++,ky+=sz,fy=py[iy]) if ((fy>=-1.0)&&(fy<=+1.0))
for (ix=0,kx=0,fx=px[ix];ix<nx;ix++,kx+=sz,fx=px[ix]) if ((fx>=-1.5)&&(fx<=+0.5))
{
// invalid qx,qy ... recompute Mandelbrot
if ((qx[ix]<0)||(qy[iy]<0))
{
for (x=0.0,y=0.0,xx=0.0,yy=0.0,i=0;(i<n)&&(xx+yy<4.0);i++)
{
q=xx-yy+fx;
y=(2.0*x*y)+fy;
x=q;
xx=x*x;
yy=y*y;
}
i=(16*i)/(n-1); if (i>16) i=16; if (i<0) i=0;
i=16-i; p[ix][iy]=i;
}
// use stored intensity
else i=p[ix][iy];
// render point with intensity i coresponding to ix,iy position in map
for (i<<=3 ,jy=0;jy<sz;jy++)
for (shd=shade8x8[i+(jy&7)],jx=0;jx<sz;jx++)
lcd.pixel(kx+jx,ky+jy,shd&(1<<(jx&7)));
}
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
lcd 和 shade8x8 内容可以在 linked SSD1306 QA 中找到。但是你可以忽略它它只是抖动和输出一个像素所以你可以直接输出 i （即使没有缩放到 <0..16>.
这里预览（在PC上因为我懒得连接相机...）：
所以它 64x32 Mandelbrot 像素显示为 128x64 抖动图像。在我的 AVR32 上，这可能比原始方法快 8 倍（可能是 3-4fps）...代码可能会更优化但是请记住 Mandelbrot 不是唯一的东西 运行 因为我有一些 ISR 处理程序处理 LCD 和我的 TTS engine based on this 的背景，从那以后我升级了很多并用于调试它（是的，它可以与渲染并行进行）。此外，我的内存不足，因为我的 3D 引擎占用了大约 11 KB 的空间（主要是深度缓冲区）。
使用此代码（内部计时器）完成预览：
static float zoom=1.0;
mandelbrot_draw<128,64,1>(+0.37,-0.1,zoom);
zoom*=1.02; if (zoom>100000) zoom=1.0;
对于非 AVR32 C++ 环境也使用这个：
//------------------------------------------------------------------------------------------
#ifndef _AVR32_compiler_h
#define _AVR32_compiler_h
//------------------------------------------------------------------------------------------
typedef int32_t S32;
typedef int16_t S16;
typedef int8_t S8;
typedef uint32_t U32;
typedef uint16_t U16;
typedef uint8_t U8;
//------------------------------------------------------------------------------------------
#endif
//------------------------------------------------------------------------------------------
[Edit2] GLSL 中更高的浮点精度
Mandelbrot 的主要问题是它需要添加指数差异很大的数字。对于 +,- 操作，我们需要对齐两个操作数的尾数，将它们添加为整数并归一化回科学记数法。但是，如果指数差异很大，则结果尾数需要的位数超过 32 位浮点数所能容纳的位数，因此仅保留 24 个最高有效位。这会产生导致像素化的舍入误差。如果您查看二进制的 32 位 float，您会看到：
float a=1000.000,b=0.00001,c=a+b;
//012345678901234567890123456789 ... just to easy count bits
a=1111101000b // a=1000
b= 0.00000000000000001010011111000101101011b // b=0.00000999999974737875
c=1111101000.00000000000000001010011111000101101011b // not rounded result
c=1111101000.00000000000000b // c=1000 rounded to 24 bits of mantissa
现在的想法是增加尾数位数。最简单的技巧是使用 2 个浮点数而不是一个：
//012345678901234567890123 ... just to easy count bits
a=1111101000b // a=1000
//012345678901234567890123 ... just to easy count b= 0.0000000000000000101001111100010110101100b // b=0.00000999999974737875
c=1111101000.00000000000000001010011111000101101011b // not rounded result
c=1111101000.00000000000000b // c=1000 rounded to 24
+ .0000000000000000101001111100010110101100b
//012345678901234567890123 ... just to easy count bits
所以结果的一部分在一个浮点数中，其余部分在另一个浮点数中...每个单个值的浮点数越多，尾数就越大。然而，在 GLSL 中将大尾数精确划分为 24 位块会很复杂且速度很慢（如果由于 GLSL 限制甚至可能的话）。相反，我们可以 select 为每个 float 指定一定范围的指数（就像上面的例子一样）。
因此在示例中，每个单个（浮点数）值有 3 个浮点数（vec3）。每个坐标代表不同的范围：
abs(x) <= 1e-5
1e-5 < abs(y) <= 1e+5
1e+5 < abs(z)
和 value = (x+y+z) 所以我们有点像 3*24 位尾数，但是范围并不完全匹配 24 位。为此，指数范围应除以：
log10(2^24)=7.2247198959355486851297334733878
而不是 10...例如这样的东西：
abs(x) <= 1e-7
1e-7 < abs(y) <= 1e+0
1e+0 < abs(z)
此外，范围 必须 selected 以便它们处理您使用的值范围 否则它将一无所获。因此，如果您的数字是 <4 范围 >10^+5 毫无意义，那么首先您需要查看您拥有的值的范围，然后将其分解为指数范围（与每个值的浮点数一样多）。
注意仍然会发生一些舍入（但比原生浮点数少得多）！！！
现在对此类数字进行运算比对普通浮点数进行运算稍微复杂一些，因为您需要将每个值作为所有分量的括号总和来处理，因此：
(x0+y0+z0) + (x1+y1+z1) = (x0+x1 + y0+y1 + z0+z1)
(x0+y0+z0) - (x1+y1+z1) = (x0-x1 + y0-y1 + z0-z1)
(x0+y0+z0) * (x1+y1+z1) = x0*(x1+y1+z1) + y0*(x1+y1+z1) + z0*(x1+y1+z1)
并且不要忘记将值规范化回定义的范围。避免添加大小 (abs) 值，因此请避免 x0+z0 等 ...
[Edit3] 新的 win32 演示 CPU 与 GPU
两个可执行文件都预设到相同的位置和缩放以显示 doubles 何时开始四舍五入。我不得不稍微升级 px,py 坐标的计算方式，因为大约 10^9 y 轴开始在这个位置偏离（对于其他位置，阈值仍然可能太大）
此处预览 CPU 与 GPU 的高倍率 (n=1600):
CPU 的 RT GIF 捕获（n=62++，GIF 缩小 4 倍）：
我已经尝试了很多算法来渲染 Mandelbrot 集，包括朴素的逃逸时间算法，以及优化的逃逸时间算法。但是，是否有更快的算法可以像我们在 YouTube 上看到的那样有效地产生真正的深度缩放。此外，我很想获得一些关于如何提高 C/C++ double
优化后的转义算法应该足够快，可以实时绘制 Mandelbrot 集。您可以使用 多线程 以便您的实施速度更快（例如，使用 OpenMP 非常容易）。如果需要，您还可以使用 SIMD 指令 手动矢量化您的代码以使其更快。如果这对您来说仍然不够快（尽管要有效地做到这一点有点复杂）。最后，您应该调整迭代次数，使其相当小。
缩放不应该对性能有任何直接影响。它只是改变了计算的输入window。但是，它确实具有 间接 影响，因为实际迭代次数会发生变化。 不应计算 window 之外的点。
双精度也应该足以正确绘制 Mandelbrot 集。但是如果你真的想要更精确的计算，你可以使用double-double precision，它提供了相当好的精度并且性能也不错。但是，手动实现 double-double precision 有点棘手，而且它仍然比仅使用 double precision 要慢得多。
我最快的解决方案通过遵循等高线边界和填充来避免在相同深度的大面积上进行迭代。有一个惩罚是可以剪掉小芽而不是绕过它们，但总而言之，为快速变焦付出的代价很小。
一个可能的效率是，如果缩放比例加倍，您已经拥有 ¼ 的点。
对于动画，我记录了每一帧的值，每次将比例加倍，并在播放时实时插入中间帧，因此动画每秒加倍一次。 double 类型允许存储超过 50 个关键帧，给出一个持续超过一分钟的动画（进入然后退出）。
实际迭代是由手工汇编程序完成的，因此一个像素完全在 FPU 中迭代。
即使是高端 CPU 与普通 GPU 相比也会慢得多。即使在 GPU 上使用简单的迭代算法，您也可以获得实时渲染。因此，在 GPU 上使用更好的算法可以获得高变焦，但是对于您需要的任何体面的算法：
- 多通道渲染，因为我们不能在 GPU 上自行修改纹理
- 高精度浮点数
floats/doubles是不够的。
这里有几个相关的问答：
- Interior Distance Estimate algorithm for the Mandelbrot set
- Perturbation theory
这可能会让你开始......
一种加快速度的方法是您可以像我在第一个 link 中那样使用 分数转义 。它提高了图像质量，同时保持较低的最大迭代次数。
第二个 link 将让您近似分形的哪些部分进出以及多远。它不是很准确，但可以用来避免计算“肯定在外面”的部分的迭代。
接下来 link 将向您展示如何获得更高的精度。
最后一个 link 是关于 扰动 这个想法是你只对一些参考点使用高精度数学，并使用它来计算它的低精度数学的邻居点在不失去精度的情况下。从未使用过，但看起来很有前途。
最后，一旦您实现了快速渲染，您可能希望以此为目标：
- How to adjust panning while zooming Mandelbrot set
这里是 GLSL 中用于单值的 3* 64 位双精度的小例子：
// high precision float (very slow)
dvec3 fnor(dvec3 a)
{
dvec3 c=a;
if (abs(c.x)>1e-5){ c.y+=c.x; c.x=0.0; }
if (abs(c.y)>1e+5){ c.z+=c.y; c.y=0.0; }
return c;
}
double fget(dvec3 a){ return a.x+a.y+a.z; }
dvec3 fset(double a){ return fnor(dvec3(a,0.0,0.0)); }
dvec3 fadd(dvec3 a,double b){ return fnor(a+fset(b)); }
dvec3 fsub(dvec3 a,double b){ return fnor(a-fset(b)); }
dvec3 fmul(dvec3 a,double b){ return fnor(a*b); }
dvec3 fadd(dvec3 a,dvec3 b){ return fnor(a+b); }
dvec3 fsub(dvec3 a,dvec3 b){ return fnor(a-b); }
dvec3 fmul(dvec3 a,dvec3 b)
{
dvec3 c;
c =fnor(a*b.x);
c+=fnor(a*b.y);
c+=fnor(a*b.z);
return fnor(c);
}
所以每个高精度值都是 dvec3 ... fnor 中的阈值可以更改为任何范围。您可以将其转换为 vec3 和 float ...
[Edit1]“快速”C++ 示例
好的，我想试试我的新
将分辨率降低 2
为抖动腾出空间，因为我的输出只是黑白
基于
使用变量最大迭代zoomn更改
n后会使最后一帧无效以强制执行完全重新计算。我知道这很慢，但在缩放范围之间的转换中只发生 3 次。最后一帧的缩放计数看起来不太好，因为它不是线性的。
可以使用最后的计数，但为此还需要记住用于迭代的复杂变量，这会占用太多内存。
记住最后一帧以及
x,y屏幕坐标映射到哪个 Mandelbrot 坐标。在每一帧计算屏幕坐标和 Mandelbrot 坐标之间的映射。
重新映射最后一帧以调整到新位置和缩放
所以只需查看来自 #3、#4 的数据，如果我们在最后一帧和实际帧中的位置相同（接近像素大小的一半），复制像素。并重新计算其余部分。
如果您的视图平滑（因此位置和缩放在每帧基础上不会有太大变化），这将极大地提高性能。
我知道它的描述有点模糊，所以这里有一个C++代码，您可以在其中推断所有疑问：
//---------------------------------------------------------------------------
//--- Fast Mandelbrot set ver: 1.000 ----------------------------------------
//---------------------------------------------------------------------------
template<int xs,int ys,int sh> void mandelbrot_draw(float mx,float my,float zoom)
{
// xs,ys - screen resolution
// sh - log2(pixel_size) ... dithering pixel size
// mx,my - Mandelbrot position (center of view) <-1.5,+0.5>,<-1.0,+1.0>
// zoom - zoom
// ----------------
// (previous/actual) frame
static U8 p[xs>>sh][ys>>sh]; // intensities (raw Mandelbrot image)
static int n0=0; // max iteraions
static float px[(xs>>sh)+1]={-1000.0}; // pixel x position in Mandlebrot
static float py[(ys>>sh)+1]; // pixel y position in Mandlebrot
// temp variables
U8 shd; // just pattern for dithering
int ix,iy,i,n,jx,jy,kx,ky,sz; // index variables
int nx=xs>>sh,ny=ys>>sh; // real Mandelbrot resolution
float fx,fy,fd; // floating Mandlebrot position and pixel step
float x,y,xx,yy,q; // Mandelbrot iteration stuff (this need to be high precision)
int qx[xs>>sh],qy[ys>>sh]; // maping of pixels between last and actual frame
float px0[xs>>sh],py0[ys>>sh]; // pixel position in Mandlebrot from last frame
// init vars
if (zoom< 10.0) n= 31;
else if (zoom< 100.0) n= 63;
else if (zoom< 1000.0) n=127;
else n=255;
sz=1<<sh;
ix=xs; if (ix>ys) ix=ys; ix/=sz;
fd=2.0/(float(ix-1)*zoom);
mx-=float(xs>>(1+sh))*fd;
my-=float(ys>>(1+sh))*fd;
// init buffers
if ((px[0]<-999.0)||(n0!=n))
{
n0=n;
for (ix=0;ix<nx;ix++) px[ix]=-999.0;
for (iy=0;iy<ny;iy++) py[iy]=-999.0;
for (ix=0;ix<nx;ix++)
for (iy=0;iy<ny;iy++)
p[ix][iy]=0;
}
// store old and compute new float positions of pixels in Mandelbrot to px[],py[],px0[],py0[]
for (fx=mx,ix=0;ix<nx;ix++,fx+=fd){ px0[ix]=px[ix]; px[ix]=fx; qx[ix]=-1; }
for (fy=my,iy=0;iy<ny;iy++,fy+=fd){ py0[iy]=py[iy]; py[iy]=fy; qy[iy]=-1; }
// match old and new x coordinates to qx[]
for (ix=0,jx=0;(ix<nx)&&(jx<nx);)
{
x=px[ix]; y=px0[jx];
xx=(x-y)/fd; if (xx<0.0) xx=-xx;
if (xx<=0.5){ qx[ix]=jx; px[ix]=y; }
if (x<y) ix++; else jx++;
}
// match old and new y coordinates to qy[]
for (ix=0,jx=0;(ix<ny)&&(jx<ny);)
{
x=py[ix]; y=py0[jx];
xx=(x-y)/fd; if (xx<0.0) xx=-xx;
if (xx<=0.5){ qy[ix]=jx; py[ix]=y; }
if (x<y) ix++; else jx++;
}
// remap p[][] by qx[]
for (ix=0,jx=nx-1;ix<nx;ix++,jx--)
{
i=qx[ix]; if ((i>=0)&&(i>=ix)) for (iy=0;iy<ny;iy++) p[ix][iy]=p[i][iy];
i=qx[jx]; if ((i>=0)&&(i<=jx)) for (iy=0;iy<ny;iy++) p[jx][iy]=p[i][iy];
}
// remap p[][] by qy[]
for (iy=0,jy=ny-1;iy<ny;iy++,jy--)
{
i=qy[iy]; if ((i>=0)&&(i>=iy)) for (ix=0;ix<nx;ix++) p[ix][iy]=p[ix][i];
i=qy[jy]; if ((i>=0)&&(i<=jy)) for (ix=0;ix<nx;ix++) p[ix][jy]=p[ix][i];
}
// Mandelbrot
for (iy=0,ky=0,fy=py[iy];iy<ny;iy++,ky+=sz,fy=py[iy]) if ((fy>=-1.0)&&(fy<=+1.0))
for (ix=0,kx=0,fx=px[ix];ix<nx;ix++,kx+=sz,fx=px[ix]) if ((fx>=-1.5)&&(fx<=+0.5))
{
// invalid qx,qy ... recompute Mandelbrot
if ((qx[ix]<0)||(qy[iy]<0))
{
for (x=0.0,y=0.0,xx=0.0,yy=0.0,i=0;(i<n)&&(xx+yy<4.0);i++)
{
q=xx-yy+fx;
y=(2.0*x*y)+fy;
x=q;
xx=x*x;
yy=y*y;
}
i=(16*i)/(n-1); if (i>16) i=16; if (i<0) i=0;
i=16-i; p[ix][iy]=i;
}
// use stored intensity
else i=p[ix][iy];
// render point with intensity i coresponding to ix,iy position in map
for (i<<=3 ,jy=0;jy<sz;jy++)
for (shd=shade8x8[i+(jy&7)],jx=0;jx<sz;jx++)
lcd.pixel(kx+jx,ky+jy,shd&(1<<(jx&7)));
}
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
lcd 和 shade8x8 内容可以在 linked SSD1306 QA 中找到。但是你可以忽略它它只是抖动和输出一个像素所以你可以直接输出 i （即使没有缩放到 <0..16>.
这里预览（在PC上因为我懒得连接相机...）：
所以它 64x32 Mandelbrot 像素显示为 128x64 抖动图像。在我的 AVR32 上，这可能比原始方法快 8 倍（可能是 3-4fps）...代码可能会更优化但是请记住 Mandelbrot 不是唯一的东西 运行 因为我有一些 ISR 处理程序处理 LCD 和我的 TTS engine based on this 的背景，从那以后我升级了很多并用于调试它（是的，它可以与渲染并行进行）。此外，我的内存不足，因为我的 3D 引擎占用了大约 11 KB 的空间（主要是深度缓冲区）。
使用此代码（内部计时器）完成预览：
static float zoom=1.0;
mandelbrot_draw<128,64,1>(+0.37,-0.1,zoom);
zoom*=1.02; if (zoom>100000) zoom=1.0;
对于非 AVR32 C++ 环境也使用这个：
//------------------------------------------------------------------------------------------
#ifndef _AVR32_compiler_h
#define _AVR32_compiler_h
//------------------------------------------------------------------------------------------
typedef int32_t S32;
typedef int16_t S16;
typedef int8_t S8;
typedef uint32_t U32;
typedef uint16_t U16;
typedef uint8_t U8;
//------------------------------------------------------------------------------------------
#endif
//------------------------------------------------------------------------------------------
[Edit2] GLSL 中更高的浮点精度
Mandelbrot 的主要问题是它需要添加指数差异很大的数字。对于 +,- 操作，我们需要对齐两个操作数的尾数，将它们添加为整数并归一化回科学记数法。但是，如果指数差异很大，则结果尾数需要的位数超过 32 位浮点数所能容纳的位数，因此仅保留 24 个最高有效位。这会产生导致像素化的舍入误差。如果您查看二进制的 32 位 float，您会看到：
float a=1000.000,b=0.00001,c=a+b;
//012345678901234567890123456789 ... just to easy count bits
a=1111101000b // a=1000
b= 0.00000000000000001010011111000101101011b // b=0.00000999999974737875
c=1111101000.00000000000000001010011111000101101011b // not rounded result
c=1111101000.00000000000000b // c=1000 rounded to 24 bits of mantissa
现在的想法是增加尾数位数。最简单的技巧是使用 2 个浮点数而不是一个：
//012345678901234567890123 ... just to easy count bits
a=1111101000b // a=1000
//012345678901234567890123 ... just to easy count b= 0.0000000000000000101001111100010110101100b // b=0.00000999999974737875
c=1111101000.00000000000000001010011111000101101011b // not rounded result
c=1111101000.00000000000000b // c=1000 rounded to 24
+ .0000000000000000101001111100010110101100b
//012345678901234567890123 ... just to easy count bits
所以结果的一部分在一个浮点数中，其余部分在另一个浮点数中...每个单个值的浮点数越多，尾数就越大。然而，在 GLSL 中将大尾数精确划分为 24 位块会很复杂且速度很慢（如果由于 GLSL 限制甚至可能的话）。相反，我们可以 select 为每个 float 指定一定范围的指数（就像上面的例子一样）。
因此在示例中，每个单个（浮点数）值有 3 个浮点数（vec3）。每个坐标代表不同的范围：
abs(x) <= 1e-5
1e-5 < abs(y) <= 1e+5
1e+5 < abs(z)
和 value = (x+y+z) 所以我们有点像 3*24 位尾数，但是范围并不完全匹配 24 位。为此，指数范围应除以：
log10(2^24)=7.2247198959355486851297334733878
而不是 10...例如这样的东西：
abs(x) <= 1e-7
1e-7 < abs(y) <= 1e+0
1e+0 < abs(z)
此外，范围 必须 selected 以便它们处理您使用的值范围 否则它将一无所获。因此，如果您的数字是 <4 范围 >10^+5 毫无意义，那么首先您需要查看您拥有的值的范围，然后将其分解为指数范围（与每个值的浮点数一样多）。
注意仍然会发生一些舍入（但比原生浮点数少得多）！！！
现在对此类数字进行运算比对普通浮点数进行运算稍微复杂一些，因为您需要将每个值作为所有分量的括号总和来处理，因此：
(x0+y0+z0) + (x1+y1+z1) = (x0+x1 + y0+y1 + z0+z1)
(x0+y0+z0) - (x1+y1+z1) = (x0-x1 + y0-y1 + z0-z1)
(x0+y0+z0) * (x1+y1+z1) = x0*(x1+y1+z1) + y0*(x1+y1+z1) + z0*(x1+y1+z1)
并且不要忘记将值规范化回定义的范围。避免添加大小 (abs) 值，因此请避免 x0+z0 等 ...
[Edit3] 新的 win32 演示 CPU 与 GPU
两个可执行文件都预设到相同的位置和缩放以显示 doubles 何时开始四舍五入。我不得不稍微升级 px,py 坐标的计算方式，因为大约 10^9 y 轴开始在这个位置偏离（对于其他位置，阈值仍然可能太大）
此处预览 CPU 与 GPU 的高倍率 (n=1600):
CPU 的 RT GIF 捕获（n=62++，GIF 缩小 4 倍）：