/*
 * Newton2.c: Newton Method 2変数のニュートン法
 */
#include <stdio.h>
#include <math.h>

double f1(double x, double y)
/*  f1(x,y) = x^2+y^2-1  */
{return x*x+y*y-1;}

double f2(double x, double y)
/*  f2(x,y) = x^3-y  */
{return x*x*x-y;}

double J11(double x, double y)
/*  df1(x,y)/dx = 2x  */
{return 2*x;}

double J12(double x, double y)
/*  df1(x,y)/dy = 2y  */
{return 2*y;}

double J21(double x, double y)
/*  df2(x,y)/dx = 3x^2  */
{return 3*x*x;}

double J22(double x, double y)
/*  df2(x,y)/dy = -1  */
{return -1;}

int ucalc(double x, double y, double *ux, double *uy)
/*
    u = -J^-1(x)f(x)
      = -1/det*( J22,-J12)*(f1(x))
               (-J21, J11) (f2(x))
*/
{
    double det;

    det = J11(x,y)*J22(x,y)-J12(x,y)*J21(x,y);
    *ux = -1/det*(J22(x,y)*f1(x,y)-J12(x,y)*f2(x,y));
    *uy =  1/det*(J21(x,y)*f1(x,y)-J11(x,y)*f2(x,y));
    return 0;
}

int Newton2(double *x, double *y)
{
    double x0, y0, ux, uy, err;
    double epsa = 1e-70, epsr = 1e-10;

    do {
        x0 = *x; y0 = *y;
        ucalc(*x,*y,&ux,&uy);
        *x = *x+ux; *y = *y+uy;
        err = hypot(*x-x0,*y-y0)/
              (epsa+epsr*(hypot(x0,y0)+hypot(*x,*y)));
        printf(" (x, y) = ( %e,  %e), err =  %e\n",*x,*y,err);
    } while (err >= 1);
    return 0;
}

/*    ------------    Main    ------------    */
int main()
{
    double x, y;
    char s[128];

    printf(" (x0, y0) = "); gets(s); sscanf(s,"%lf %lf",&x,&y);
    Newton2(&x,&y);
    putchar('\n');
    printf(" (x, y) = ( %e,  %e)\n", x, y);
    return 0;
}
/*   ---   実行例: 新濃・船田「数値解析の基礎」 p.56 問1   ---   */
/*
Z:\src\C>Newton2
 (x0, y0) = 1 0.5
 (x, y) = ( 8.500000e-01,  5.500000e-01), err =  7.421595e+08
 (x, y) = ( 8.266083e-01,  5.634235e-01), err =  1.339916e+08
 (x, y) = ( 8.260317e-01,  5.636240e-01), err =  3.051883e+06
 (x, y) = ( 8.260314e-01,  5.636242e-01), err =  1.809697e+03
 (x, y) = ( 8.260314e-01,  5.636242e-01), err =  5.895258e-04

 (x, y) = ( 8.260314e-01,  5.636242e-01)
*/