0

今晚我很想回家，使用 Mandelbrot/Julia 分形发生器。这是我之前在 C++ 中进行的一个项目，但这次我将在 C# 中试一试，这样我就可以更轻松地将其制作为 WinForms 应用程序，使用多线程生成图像等。

然而，我考虑使用但看起来可能最终会成为拐杖的一件事是内置的 Complex 结构。乍一看，看起来很棒；我需要的所有操作都内置到一个结构中。但是，仔细观察会发现该结构使用双精度数。双打只有 15-16 sig figs；他们的力量是他们的范围。我只对复平面上 (-2,-2) 到 (2,2) 的范围感兴趣，所以我对范围不感兴趣；精度。小数有 28 位有效数字。

因此问题；去这里的路是什么？

  1. 使用内置的 Complex 类型；毕竟这只是一个附带项目。
  2. 使用小数而不是双精度数滚动我自己的 DecimalComplex 类型。似乎是一条不错的“中间道路”；我需要的数学并不难（乘法和加法），它几乎可以使可用精度翻倍。
  3. 孤注一掷; 实现一个 BigDecimal，利用 BigInteger 并自己跟踪小数位，并将其用于 UberComplex 结构。IIRC 上次我搞砸这个时，我使用了类似的想法，但使用了无符号的 64 位长，这让我达到了一个缩放级别，如果整个集合以 360dpi 绘制，它需要一个大致区域的表面的太阳系。我可能会做得比这更好。
4

1 回答 1

0

如果这是您第一次编写 Mandelbrot 生成器，我建议您使用该Complex类型。过早的优化是万恶之源；而且，Mandelbrot 的一半酷之处在于您可以在几行代码中生成令人惊叹的图像。花两个小时来实现你自己的BigDecimal那种失败的目的。

但是，我认为您没有提到另一种选择：使用定点数学代替 IEEE 浮点数。定点运算对于 Mandelbrot 生成确实很方便，因为加法完全一样，乘法几乎一样：(x * y)变成(x * y) >> MANTISSABITS.

剧透警告！这是我几年前写的一个定点 Mandelbrot 生成器。它使用小数点后的 28 位。如果您的工具链支持int128_t，或者您不介意进入汇编程序，或者您不介意长时间进行乘法运算，则可以轻松替换int32_t为int64_t并获得 60 位分辨率。

#include <stdio.h>
#include <stdint.h>
#include "ImageFmtc.h"

/* Represent each component of the complex number
 * as a fixed-point value between -4 and +4. There
 * are 4 bits before the binary point and 28 bits after. */
typedef int32_t REALFIXED;
typedef int32_t IMAGFIXED;

#define MandelWidth   300
#define MandelHeight  300

unsigned char Image[MandelWidth][MandelHeight];


void mandelbrot(REALFIXED centerx, IMAGFIXED centery, REALFIXED windowsize)
{
    int x, y;
    unsigned char *image_ptr = &Image[0][0];

    for (y = 0; y < MandelHeight; ++y) {

    const IMAGFIXED z_imag = (1LL * windowsize * y / MandelHeight) + centery;

    for (x = 0; x < MandelWidth; ++x) {

        const REALFIXED z_real = (1LL * windowsize * x / MandelWidth) + centerx;

        REALFIXED t_real = z_real;
        IMAGFIXED t_imag = z_imag;

        int iter;

        for (iter=0; iter < 255; ++iter) {
        long long r2 = (1uLL * t_real * t_real) >> 28;
        long long i2 = (1uLL * t_imag * t_imag) >> 28;
        long long two_ri = (1uLL * t_real * t_imag) >> 27;
        if (r2 + i2 > (4LL << 28))
            break;
        t_real = (r2 - i2) + z_real;
        t_imag = (two_ri) + z_imag;
        }

        *image_ptr++ = iter & 0xFF;
    }
    }
    return;
}

REALFIXED double2fixed(double x)
{
    return (REALFIXED)(x * (1uLL << 28));
}

int main(int argc, char **argv)
{
    double cx, cy, scale;
    sscanf(argv[1], "%lf", &cx);
    sscanf(argv[2], "%lf", &cy);
    sscanf(argv[3], "%lf", &scale);
    mandelbrot(double2fixed(cx-scale/2), double2fixed(cy-scale/2), double2fixed(scale));
    WritePGM5("mandel.pgm", (unsigned char *)Image, MandelWidth, MandelHeight);
    return 0;
}

关于太阳系的大小：即使您的数字可能具有 60 位的分辨率，但这并不意味着您实际上可以计算出以 60 位分辨率设置的 Mandelbrot 像素精确渲染。乘法上的舍入误差意味着当你放大时，你的图像会变得比你预期的要快得多。我确定有人做过数学计算，但我不知道它是不是在我的脑海中。关键字“错误传播”。

于 2012-09-20T23:12:02.083 回答