//     ニュートン法による非線形方程式の解法   

#include<stdio.h>

#include<math.h>

#define EPSILON 1.0E-7         //  反復におけるｘの変動幅の縮小限界 

double a[20];                  //　係数                  

double r,root;                 //  r:近似値初期値　      

int    n,k;                    //  n次方程式 ,k:反復回数 

double f(double),df(double);   //  df:df/dx              

int    sw;                     //  sw=0 : no convergent  

void main()

{

    char ym,yn;

    int i;

    double Newton(void);

    printf("\n\n------------------Newton Method------------------\n\n");

    printf("方程式の次数 : n = ?");

    scanf("%d",&n);

    printf("係数(高次から) : a[i] = ? (i=0,,,%2d)\n",n);

        printf("各値は、空欄で区切る");

    for(i=0;i<=n;i++)

                scanf("%lf",&a[i]);

www:sw=1;

    printf("解の近似計算の初期値 : r = ?");

        scanf("%lf",&r);

        printf("\n");

    if(f(r)==0)

                printf(" root = %f\n\n",r);

    else

    {

                if(df(r)<0.0001)

                {

                        printf("df:%lf",df(r));

                        printf("初期値ｒが不適当です.\n");

                        printf("No convergent. because f'(r)=0 or f'(r)is very small.\n\n");

                }

                else

                {

                        if(sw)

                        {

                                root = Newton();

                                printf("公差＜%10.8lfで反復計算を打ち切りました.\n\n",EPSILON);

                                printf("収束値（解）: root = %f\n\n",root);

                        }

                        else

                                printf("no convergent !\n\n");

                }

    }

    printf("別の初期値で再度計算しますか.<y/n>?");

    ym=getchar();

        yn=getchar();

    if(yn=='y'||yn=='Y') goto www;

}

double Newton()

{

    double d,dd,x;

    k=0;

    do

    {

                k++;

                x=r-f(r)/df(r);

                d=fabs(x-r);

                printf("k=%2d : x=%11.8lf\n",k,x);

                r=x;

                if(k>1)

                    if(d>dd)

                        {

                                sw=0;

                                break;

                        }

                        dd=d;

        }

    while(d>EPSILON);

    return x;

}

double f(double x)

{

    double b[20];

    int i;

    b[0]=a[0];

    for(i=1;i<=n;i++)

                b[i]=b[i-1]*x+a[i];

    return b[n];

}

double df(double x)

{

    double c[20];

    int i;

    c[0]=n*a[0];

    for(i=1;i<n;i++)

                c[i]=c[i-1]*x+(n-i)*a[i];

    return c[n-1];

}