可变指数类Mandelbrot集C++计算优化技术咨询

我正在C++中优化广义Mandelbrot集的计算，核心差异是指数不是固定的2.0，而是从输入读取的可变值（如2.33、2.11351等），迭代逻辑等效于：

z = pow(z, powerReal) + c;

仅需统计未逃逸点的数量，无需生成图像。

当前编译命令：

g++ --std=c++17 -Wall -Wextra -fopenmp -march=native -O3

perf分析显示耗时主要集中在sincos、atan2、cexp、clog、hypot、log、exp等libm函数上，希望针对以下方向优化：

1. 是否存在更快的复数z^powerReal计算方式？
2. 能否在保持数值行为一致的前提下，规避部分复数/对数/指数运算开销？
3. 此场景下SIMD是否可行？非整数指数是否会造成阻碍？
4. 这类广义Mandelbrot集是否有安全的数学捷径？

当前计算耗时约0.7s，目标是优化到0.15s左右，且必须保证numberInSet结果完全一致，是否可通过预计算等方式优化？

附当前实现代码：

#include <complex>
#include <cstdio>
#include <iostream>

constexpr int max_iter = 223;

constexpr int x_resolution = 1923;
constexpr int y_resolution = 1447;

constexpr double x_step = 4.0 / x_resolution;
constexpr double y_step = 4.0 / y_resolution;

constexpr int x_start = 100;
constexpr int x_end = 1305;
constexpr int y_start = 230;
constexpr int y_end = 1417;

int mandelbrot_draw(double powerReal) {
    int numberInSet = 0;

#pragma omp parallel for collapse(2) reduction(+ : numberInSet) schedule(dynamic, 16)
    for (int py = y_start; py < y_end; py++) {
        for (int px = x_start; px < x_end; px++) {
            const double y0 = 2.0 - py * y_step;
            const double x0 = -2.0 + px * x_step;

            std::complex<double> z(0.0, 0.0);

            unsigned short iteration = 0;

            for (; iteration < max_iter; iteration++) {
                const std::complex<double> logZ = std::log(z);

                const std::complex<double> scaledLogZ(
                    powerReal * logZ.real(),
                    powerReal * logZ.imag()
                );

                z = std::exp(scaledLogZ);

                z.real(z.real() + x0);
                z.imag(z.imag() + y0);

                if (z.real() * z.real() + z.imag() * z.imag() >= 4.0) {
                    ++iteration;
                    break;
                }
            }

            numberInSet += (iteration == max_iter);
        }
    }

    return numberInSet;
}

int main() {
    double powerReal;

    // The input is a seed fraction, for example:
    // input "33" means powerReal = 2.33
    std::string seedFraction;
    std::cin >> seedFraction;

    int value = 0;
    double divisor = 1.0;

    for (char c : seedFraction) {
        value = value * 10 + (c - '0');
        divisor *= 10.0;
    }

    powerReal = 2.0 + value / divisor;

    int numberInSet = mandelbrot_draw(powerReal);

    std::printf("%d\n", numberInSet);
}

1. 更快的复数z^powerReal计算方式 #

当前代码用log(z) * powerReal再exp的逻辑，本质是极坐标转换：z^k = r^k * e^(i*k*theta)（r = |z|，theta = arg(z)）。直接手动拆分极坐标计算比调用std::log+std::exp更高效，因为标准库函数会处理分支切割等额外逻辑，而我们可以聚焦核心计算：

// 替换原log+exp逻辑
double r = std::hypot(z.real(), z.imag());
double theta = std::atan2(z.imag(), z.real());
z = std::polar(std::pow(r, powerReal), powerReal * theta);

std::polar是标准库专门用于从模长和辐角构造复数的函数，内部实现比手动计算sin/cos更高效，同时能保证和原代码的数值一致性。

2. 保持数值一致前提下规避运算开销 #

跳过冗余初始迭代 #

初始z为(0,0)，第一次迭代的pow(0, powerReal)结果还是0，加c后等于c，完全可以跳过这次无意义的循环：

unsigned short iteration = 1;
std::complex<double> z(x0, y0);
double z_mag_sq = x0*x0 + y0*y0;

for (; iteration < max_iter; iteration++) {
    // 计算z^powerReal...
    z += std::complex<double>(x0, y0);
    z_mag_sq = std::norm(z);
    if (z_mag_sq >= 4.0) {
        ++iteration;
        break;
    }
}

复用模长平方计算 #

逃逸判断需要的|z|²可以在每次迭代后保存，避免重复计算两次乘法和一次加法，用std::norm(z)也能简化代码。

预计算所有c值 #

将所有像素对应的c = (x0, y0)提前计算并存入数组，避免嵌套循环中重复计算坐标转换：

std::vector<std::complex<double>> cs;
cs.reserve((y_end - y_start) * (x_end - x_start));
for (int py = y_start; py < y_end; py++) {
    double y0 = 2.0 - py * y_step;
    for (int px = x_start; px < x_end; px++) {
        double x0 = -2.0 + px * x_step;
        cs.emplace_back(x0, y0);
    }
}

// 后续并行遍历cs数组即可

3. SIMD可行性分析 #

非整数指数确实增加了SIMD的实现难度，但并非不可行：

• 核心思路：利用SIMD寄存器批量处理4个（AVX2）或8个（AVX-512）独立的点，同时计算它们的模长、辐角、幂次和逃逸判断。
• 依赖指令：使用AVX系列的SIMD数学函数，比如_mm256_pow_pd、_mm256_atan2_pd、_mm256_cos_pd、_mm256_sin_pd，这些指令的精度和标量libm函数一致，能保证结果不变。
• 注意事项：需要用掩码标记活跃元素（未逃逸的点），避免对已逃逸的点做无效计算。手动编写SIMD代码复杂度较高，但能带来3-4倍的性能提升，是达到0.15s目标的关键。

此外，可以添加编译选项-fopenmp-simd，让编译器尝试自动对循环做SIMD优化，降低手动实现的成本。

4. 广义Mandelbrot集的数学捷径 #

对于非整数指数的广义Mandelbrot集，没有像标准集那样的通用逃逸公式，但可以做以下安全优化：

• 预计算常量因子：powerReal是单次运行的固定值，无需在迭代循环中重复传递或计算相关常量。
• 严格边界判断：仅使用|z|² >=4作为逃逸条件，不要使用近似的提前逃逸判断（比如|z| > |c|+1），避免破坏数值一致性。
• 分支切割一致性：保持atan2的辐角范围（[-π, π]）与原代码std::log的分支逻辑一致，确保幂次计算的结果完全匹配。

内容的提问来源于stack exchange，提问作者Timon Kobusch