python - 复数和分形

我正在尝试将此处的代码http://code.activestate.com/recipes/577166-newton-fractals/改编为C，并且遇到了一些麻烦。我正在使用C99的复杂类型。

我基本上已经尝试了几种无效的方法。在每个像素上，每次都会一直进行到最大迭代，因此图像显示为纯色。

关于C中类型的工作方式，我是否有根本上的错误，似乎应该起作用，我完全准确地重构了算法。

    //newt.c

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

complex f(complex z);


const int WIDTH = 512, HEIGHT = 512;

const double SCALED_X_MAX = 1.0;
const double SCALED_X_MIN = -1.0;
const double SCALED_Y_MAX = 1.0;
const double SCALED_Y_MIN = -1.0;



const int MAX_ITERATIONS = 20;
const int EPSILON = 1e-3;

int main(int argc, char **argv) {
    const double SCALED_WIDTH = SCALED_X_MAX - SCALED_X_MIN ;
    const double SCALED_HEIGHT = SCALED_Y_MAX - SCALED_Y_MIN ;

    FILE * image = fopen("newton.ppm", "w") ;
    fprintf(image, "P3\n");
    fprintf(image, "%d %d\n" , WIDTH , HEIGHT ) ;
    fprintf(image, "%d\n", 255) ;
    for ( int y = 0 ; y < HEIGHT ; ++y) {
        double zy = y * (SCALED_HEIGHT)/(HEIGHT-1) + SCALED_Y_MIN;
        for ( int x = 0 ; x < WIDTH ; ++x ) {
            double zx = x * (SCALED_WIDTH)/(WIDTH-1) + SCALED_X_MIN;
            complex z = zx + zy*I;
            int iteration = 0;
            while(iteration < MAX_ITERATIONS ) {
//                complex h=sqrt(DBL_EPSILON) + sqrt(DBL_EPSILON)*I;
                double h=1e-6;
/*
                complex zph = z + h;
                complex dz = zph - z;
                complex slope = (f(zph) - f(z))/dz;
*/
                complex volatile dz = (f(z + (h+h*I))  - f(z)) / (h+I*h) ;
                complex z0 = z - f(z) / dz;
                //fprintf(stdout,"%f \n", cabs(z0 - z ));
                if ( cabs(z0 - z) < EPSILON){
                    break;
                }
                z = z0;
                iteration++;
            }
            if (iteration != MAX_ITERATIONS) fprintf(stdout, "%d " , iteration );
            fprintf(image,"%3d %3d %3d ", iteration % 4 * 64 ,
                                       iteration % 8 * 32 ,
                                       iteration % 16 * 16);
        }
       fprintf(image, "\n") ;
    }
    fclose(image) ;

    exit(0);

}
complex f(complex z ) {
    return cpow(z,3)-1.0 ;
}

最佳答案

在检查了复杂的数学并且没有看到任何问题之后，我注意到您的错误仅是由于该行中的整数除法引起的

zy = y * (SCALED_HEIGHT)/(HEIGHT-1) + SCALED_Y_MIN;

对于zx同样如此。因为（SCALED_HEIGHT）和（HEIGHT-1）都是整数，所以不会为您提供所需的浮点结果。尝试使用：

zy = y * SCALED_HEIGHT * 1.0/(HEIGHT-1) + SCALED_Y_MIN;

对于zx同样如此。

编辑：对不起，上面是错误的。您的SCALED_HEIGHT实际上是原来的两倍，因此上述情况实际上还可以。真正的问题很简单

const int EPSILON = 1e-3;

实际上，这总是返回零，因为它是整数。您需要使EPSILON为浮点类型。