我的Mandelbrot集合代码哪里出了问题？

Question

我正在尝试用C语言实现Mandelbrot集合，但是遇到了一个奇怪的问题。我的代码如下：

#include <stdio.h>
#include <math.h>
#include <complex.h>

int iterate_pt(complex c);

int main() {
FILE *fp;
fp = fopen("mand.ppm", "w+");


double crmin = -.75;
double crmax = -.74;
double cimin = -.138;
double cimax = -.75; //Changing this value to -.127 fixed my problem.

int ncols = 256;
int nrows = 256;
int mand[ncols][nrows];
int x, y, color;
double complex c;

double dx = (crmax-crmin)/ncols;
double dy = (cimax-cimin)/nrows;

for (x = 0; x < ncols; x++){
    for (y = 0; y < nrows; y++){
        double complex imaginary = 0+1.0i;
        c = crmin+(x*dx) + (cimin+(y*dy)) * imaginary;
        mand[x][y] = iterate_pt(c);
    }
}

printf("Printing ppm header.");
fprintf(fp, "P3\n");
fprintf(fp, "%d %d\n255\n\n", ncols, nrows);

for (x = 0; x < ncols; x++) {
    for (y = 0; y < nrows; y++){
        color = mand[x][y];
        fprintf(fp, "%d\n", color);
        fprintf(fp, "%d\n", color);
        fprintf(fp, "%d\n\n", color); //Extra new line added, telling the ppm to go to next pixel.
    }
}
fclose(fp);

return 0;
}

int iterate_pt(double complex c){
double complex z = 0+0.0i;
int iterations = 0;
int k;
for (k = 1; k <= 255; k++) {
    z = z*z + c;
    if (sqrt( z*conj(z) ) > 50){
        break;
    }
    else
        ++iterations;
}
return iterations;
}

然而，这个程序的输出结果是以ppm文件存储的，看起来像这样：。感谢您的帮助！

- Nathan Jones

---

请根据此伪代码检查您的逻辑：http://en.wikipedia.org/wiki/Mandelbrot_set#For_programmers - Blender

我正在尝试用不同的技术来表示实数和虚数，并忘记移除它了。 - Nathan Jones

我不明白你为什么这样做：if (sqrt( z*conj(z) ) > 50){。如果你想要获取幅值，可以使用 cabs(z)。 - Tom Zych

发帖提出一个带有几何内容的编程问题，获得额外加分 +1。 - Peter

关于我之前的评论：傻了吧唧的。sqrt(z * conj(z)) 就是幅度，如果我花十秒钟想一想就会意识到这一点。真正的改进应该是跳过sqrt并将其与2500进行比较。 - Tom Zych

Answer 1

尝试将cimax设置为-0.127，我也在这个项目上工作，看起来它似乎起到了作用 ;)

- geekingreen

Answer 2

代码看起来不错。但是你的起始矩形不正确！ 你正在使用：

Real ranage [  -.75 ,  -.74 ]
Imag range  [ -.138 ,  -.75 ]

您确定这是您想要的吗？对我来说，这个y轴伸展得太过夸张了。

另外，标准的曼德博集算法通常使用

magnitude > 2

而不是50。 作为一个逃生检查。虽然这不应该影响集合的实际形状。

- Sanjay Manohar

---

搞定了！将比例尺更改为类似于-.127的值给了我想要的结果。哦，那是一个非常愚蠢的错误。 - Nathan Jones

Answer 3

顺便说一下，计算z*conj(z)的平方根是没有意义的。只需将不等式两边都平方，得到if (z*conj(z) > 2500)，这样可以提高性能。

- Jens