2

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 ;
}
4

1 に答える 1

3

複雑な数学をチェックして問題が見られなかった後、あなたのエラーは単に行の整数除算が原因であることに気づきました

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は実際には2倍だったので、上記は実際には問題ありませんでした。本当の問題は単に列に並んでいる

const int EPSILON = 1e-3;

これは整数であるため、実際には常にゼロを返します。EPSILONを浮動小数点型にする必要があります。

于 2013-02-25T09:33:33.917 に答える