/*
    newton2_4.c: 別の実行例
    x^2 + (y/2)^2 = 1, y^2 = x^2 + 1 -> (x,y)=(0.7746,1.265)
 */
#include <stdio.h>
#include <math.h>
#define ND 3

int func2(double x[], double f[])
{
    f[1] = x[1] * x[1] + x[2] * x[2] / 4 - 1; //  f1(x,y) = x^2 + (y/2)^2 - 1
    f[2] = x[1] * x[1] - x[2] * x[2] + 1;     //  f2(x,y) = x^2 - y^2 + 1
    return 0;
}

int Jcalc2(double x[], double J[][ND])
{
    J[1][1] = 2 * x[1];  //  df1(x,y)/dx = 2x
    J[1][2] = x[2] / 2;  //  df1(x,y)/dy = y/2
    J[2][1] = 2 * x[1];  //  df2(x,y)/dx = 2x
    J[2][2] = -2 * x[2]; //  df2(x,y)/dy = -2y
    return 0;
}

int func(double x[], double f[])
{
    f[1] = x[1] * x[1] + x[2] * x[2] - 1; //  f1(x,y) = x^2 + y^2 - 1
    f[2] = x[1] * x[1] * x[1] - x[2];     //  f2(x,y) = x^3 - y
    return 0;
}

int Jcalc(double x[], double J[][ND])
{
    J[1][1] = 2 * x[1];        //  df1(x,y)/dx = 2x
    J[1][2] = 2 * x[2];        //  df1(x,y)/dy = 2y
    J[2][1] = 3 * x[1] * x[1]; //  df2(x,y)/dx = 3x^2
    J[2][2] = -1;              //  df2(x,y)/dy = -1
    return 0;
}

int ucalc(double x[], double u[], int (*func)(double*, double*),
 int (*Jcalc)(double*,double J[][ND]))
/*
    u = -J^-1(x)*f(x)
 */
{
    double det, f[ND], J[ND][ND];

    func(x, f);
    Jcalc(x, J);
    det = J[1][1] * J[2][2] - J[1][2] * J[2][1];
    u[1] = -1 / det * (J[2][2] * f[1] - J[1][2] * f[2]);
    u[2] =  1 / det * (J[2][1] * f[1] - J[1][1] * f[2]);
    return 0;
}

int newton2(double x[], int (*func)(double*, double*),
 int (*Jcalc)(double*,double J[][ND]), double eps)
{
    int n = 0;
    double u[ND], err;

    printf("# n, x, y, err\n");
    printf("%4d, % .15e, % .15e\n", n, x[1], x[2]);
    do {
        n++;
        ucalc(x, u, func, Jcalc);
        x[1] += u[1]; x[2] += u[2];
        err = hypot(u[1],u[2]);
        printf("%4d, % .15e, % .15e, % .15e\n", n, x[1], x[2], err);
    } while (err >= eps);
    return 0;
}

//  ------------    Main    ------------
int main(void)
{
    double x[ND], eps = 1e-10;
    char s[128];

    x[1] = 1; x[2] = 0.5;
    fprintf(stderr,"# x0 y0 = "); fgets(s,128,stdin);
    sscanf(s,"%lf%lf",&x[1],&x[2]);
    newton2(x, func2, Jcalc2, eps);
    printf("\n# (x, y) = (% .15e, % .15e)\n", x[1], x[2]);
    return 0;
}
/*  ---   実行例   ---
Z:\proen1\kadai4>newton2_4
# x0 y0 = 1 0.5
# n, x, y, err
   0,  1.000000000000000e+000,  5.000000000000000e-001
   1,  8.000000000000000e-001,  1.850000000000000e+000,  1.364734406395618e+000
   2,  7.750000000000000e-001,  1.357432432432433e+000,  4.932015902442228e-001
   3,  7.745967741935483e-001,  1.268064150511886e+000,  8.936919158457984e-002
   4,  7.745966692414905e-001,  1.264914984197939e+000,  3.149166315695276e-003
   5,  7.745966692414833e-001,  1.264911064073426e+000,  3.920124513106996e-006
   6,  7.745966692414833e-001,  1.264911064067352e+000,  6.074514284039246e-012

# (x, y) = ( 7.745966692414833e-001,  1.264911064067352e+000)
*/