Mandelbrot分形

Question

这是一个输出具有mandelbrot分形的.ppm文件的代码。

我怎么才能优化这个？

#include<bits/stdc++.h>
using namespace std;

int findMandelbrot(double cr, double ci, int max_iterations)
{
    int i = 0;
    double zr = 0.0, zi = 0.0;
    while (i < max_iterations && zr * zr + zi * zi < 4.0)
    {
        double temp = zr * zr - zi * zi + cr;
        zi = 2.0 * zr * zi + ci;
        zr = temp;
        ++i;
    }
    return i;
}

double mapToReal(int x, int imageWidth, double minR, double maxR)
{
    double range = maxR - minR;
    return x * (range / imageWidth) + minR;
}

double mapToImaginary(int y, int imageHeight, double minI, double maxI)
{
    double range = maxI - minI;
    return y * (range / imageHeight) + minI;
}

int main()
{
    ifstream f("input.txt");
    int imageWidth, imageHeight, maxN;
    double minR, maxR, minI, maxI;

    if (!f)
    {
        cout << "Could not open file!" << endl;
        return 1;
    }

    f >> imageWidth >> imageHeight >> maxN;
    f >> minR >> maxR >> minI >> maxI;

    ofstream g("output_image.ppm");
    g << "P3" << endl;
    g << imageWidth << " " << imageHeight << endl;
    g << "255" << endl;


    double start = clock();

    for (int i = 0; i < imageHeight; i++)
    {
        for (int j = 0; j < imageWidth; j++)
        {
            double cr = mapToReal(j, imageWidth, minR, maxR);
            double ci = mapToImaginary(i, imageHeight, minI, maxI);

            int n = findMandelbrot(cr, ci, maxN);

            int r = ((int)sqrt(n) % 256);
            int gr = (2*n % 256);
            int b = (n % 256);

            g << r << " " << gr << " " << b << " ";
        }
        g << endl;

        if(i == imageHeight / 2) break;
    }

    cout << "Finished!" << endl;

    double stop = clock();

    cout << (stop-start)/CLOCKS_PER_SEC;
    return 0;
}

我一直走到imageHeight /2，因为在photoshop里，我可以复制另一半。

我在考虑逻辑的力量，但我尝试了一些东西，并且只使用整数.最后，我展示了运行一些数据所需的时间。我尝试了一些互联网上的东西，但什么也没成功，也许我没有把它纠正过来。

这是具有以下输入的输出: 512 512 512 -1.5 0.7 -1.0

花了0.315秒才生产出这个。

Answer 1

您的代码生成一个“普通”的PPM文件，其中所有示例值都表示为ASCII十进制数(用空格分隔)。因此，大量的时间用于格式化和写入数据。

另一种方法是编写“二进制”PPM文件。它有神奇的"P6“而不是"P3"，所有示例值都表示为单个字节。(详见PPM格式规范。)

因此，您可以创建文件并使用(使用更多描述性变量名称的文件流)编写标题：

ofstream ppmFile("output_image.ppm", ios::out | ios::binary);
ppmFile << "P6" << endl;
ppmFile << imageWidth << " " << imageHeight << endl;
ppmFile << "255" << endl;

(endl除了编写换行符外，还刷新输出文件，这在这里是不必要的。另一方面，这并不有害，因为这个部分并不是性能关键。)

在内环中，一个RGB三重奏被附加到

int n = findMandelbrot(cr, ci, maxN);

uint8_t rgbTriplet[3];
rgbTriplet[0] = (int)sqrt(n);
rgbTriplet[1] = 2 * n;
rgbTriplet[2] = n;
ppmFile.write((char *)rgbTriplet, 3);

在每一行之后都不会写换行符。将int截断为无符号类型uint8_t是很好的定义，因此不需要显式余数操作% 256。

在我对MacBook的测试中，这减少了输入数据"512、512、512 -1.5、0.7 - 1.0“的运行时间，从≈0.25秒减少到≈0.14秒。

Answer 2

我同意@user140417所说的大部分内容(尽管我同意对pow的调用时使用Deduplicator )。因此，严格遵守性能，我有以下建议。

预先计算所有你能做的事情--

在main()中，您正在计算内环的每一次迭代中的实部和虚部。这太昂贵了。每一行都有相同的实数，每一列都有相同的假想值。因此，在开始主回路之前，先预先计算1行雷亚尔和1行想象力。

你的肤色也是一样。因为您知道maxN，所以可以有3个数组--一个用于红色通道，一个用于绿色通道，另一个用于蓝色通道。只需预先计算所有的maxN可能值。然后在内部循环中，它只是一个表查找，而不是一个计算。

使用SIMD技术

你的数学大部分是为不同的浮点数字计算相同的循环。这正是SSE、AVX等SIMD指令集发明的目的。您可以计算每条指令的2或4个值，从而大大提高实现的速度。你还可以用…做进一步的研究

多线程

如果将输入区域分解为块或条，则可以有多个线程，每个线程同时在不同的块或条上工作。在多核机器上，这将以近乎线性的方式(根据我的经验)加速到一定程度。把它和SIMD指令结合起来，你就能真正地让东西动起来了！

使用GPU

如果您真的想快速完成任务，可以使用OpenGL计算着色器、OpenCL内核或其他您喜欢使用的库编写GPU实现。GPU通常有数百到数千个浮点处理器的核心，这使得它非常适合。

Answer 3

可读性

不包括bits/stdc++.h

这是一个不可移植的头文件，超出了需要，增加了编译时间和二进制文件的大小。

不使用using namespace std;

这与上述情况结合在一起，极大地污染了全局命名空间。如果您决定使用第三方库(用于精确的数学、图像绘制或其他什么)，则很有可能会遇到冲突。

使用更多的描述性变量名称

所有的zr，zi，等。不仅很难阅读，而且很容易做一个很难调试的致命错误。您还可以使用std::complex，或者至少使用std::pair。如果以后决定添加缩放或平移功能，如果继续使用这些类型的变量名，代码的可读性将迅速降低。

改名mapToReal和mapToImaginary

也许你使用的是数学上正确的定义，但从我所知道的情况来看，这个函数的更好的名称是“规范化”或“比例”，因为“地图”太抽象/模糊。

Optimization

因子在findMandelbrot

中的算法

您总共要执行(我数)6次乘法，如果稍微重构代码，则可以减少到3次：

std::complex z{0, 0};
std::complex square_temp{0, 0};
while (square_tmp.real + square_tmp.imag <= 4.0) {
    z.imag = std::pow(z.real + z.imag, 2.0) - square_temp.real - square_temp.imag;
    z.imag += c.imag;
    z.real = square_temp.real - square_temp.imag + c.real;
    square_temp.real = std::pow(z.real, 2);
    square_temp.imag = std::pow(z.imag, 2);
}
// Note code untested

考虑使用线程/工作人员

如果将图像块传递给多个线程，则有可能获得速度增长。你是如何决定如何实现这一点的。