実際のプログラム

以下にそのプログラムを挙げておきます. プログラムが長いので, 2ペイジ わたっております. ニュートン法の漸化式の実部, 虚部は手計算です. プログラムの中で, 函数 hantei は, 収束までの回数を返しています.

/* File name 8-1.c */
#include <stdio.h>
#include <math.h>

#define BASE 255
#define SIZE 512
#define BOUNDS 1.0
#define EPSILON 1.0E-7
#define DUMMY 2.0

int hantei(double x, double y);
double cabs(double x, double y);

main()
{
    double interval=2*BOUNDS/SIZE;
    double x; 
    double y= BOUNDS;
    int   i, j, result;

    printf("P3\n%d %d\n%d\n", SIZE, SIZE, BASE);

    for(i=0; i < SIZE; i++){
        x = -BOUNDS;
        for(j=0; j < SIZE; j++){
            result=hantei(x,y);
            if (result==1)
                printf("%d %d %d\n", 0, 0, BASE);
            else if (result == 2)
                printf("%d %d %d\n", 0, BASE, 0);
            else if (result == 3)
                printf("%d %d %d\n", BASE, 0, 0);
            else 
                printf("%d %d %d\n", 0, 0, 0); 
            x += interval;
        }
        y -= interval;
    }
}

int hantei(double x, double y)
{
    double u=DUMMY;
    double v=DUMMY;
    double nextu=x;
    double nextv=y;
    double denom;

    if (x==0 && y==0) return 0;

    while(cabs(nextu-u, nextv-v)>EPSILON){
        u=nextu;
        v=nextv;
        denom=3*(u*u*u*u+2*u*u*v*v+v*v*v*v);
        nextu=(2.0*u/3.0) + (u*u-v*v)/denom;
        nextv=(2.0*v/3.0)-2*u*v/denom;
    }
    if (fabs(u-1.0)< 0.5)
        return 1;
    else if (fabs(v-(sqrt(3)/2.0))<0.5)
        return 2;
    else if (fabs(v+(sqrt(3)/2.0))<0.5)
        return 3;
    else
        return 0;
}

double cabs(double a, double b)
{
    return sqrt(a*a+b*b);
}

プログラムが完成したら, 実行して見て下さい. 多量の出力が出ますので, more にリダイレクションします. 次の様になっていれば間違い無いでしょう.

e023199@cc> ./a.out|more
P3
512 512
255
0 255 0
0 255 0
0 255 0
0 255 0
0 255 0
.........

プログラムは 標準出力にデータを出力しますから, ファイルに保存するには リダイレクションを使います.

e023199@cc> ./a.out > cnewton.ppm

図形を見るには, ASTEC-X で login して, ImageView でみます. root ウィンドウで右ボタンを押して, アプリケーションから ImageView を選びます. アプリケーションが起動したらファイルメニューから, cnewton.ppm というファイルを選択します.