/*
    newton2_1.c: f(x), J(x) をまとめたバージョン
    x^2 + y^2 = 1, y = x^3 -> (x, y) = (0.826, 0.5636), (-0.826, -0.5636)
 */
#include <stdio.h> // printf, fprintf, fgets, sscanf
#include <math.h>  // hypot(hypotenuse): 直角三角形の斜辺を求める関数
#define ND 3

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

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

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

    func(x, y, f);  // 配列f[]に f(x,y) の各成分を格納
    Jcalc(x, y, J); // 配列J[][]に J(x,y) の各成分を格納
    det = J[1][1] * J[2][2] - J[1][2] * J[2][1];
    *ux = -1 / det * (J[2][2] * f[1] - J[1][2] * f[2]);
    *uy =  1 / det * (J[2][1] * f[1] - J[1][1] * f[2]);
    return 0;
}

int newton2(double *x, double *y, double eps)
{
    int n=0;
    double ux, uy, err;

    printf("# n, x, y, err\n");
    printf("%4d, % .15e, % .15e\n", n, *x, *y);
    do {
        n++;
        ucalc(*x, *y, &ux, &uy);
        *x += ux; *y += uy; // x_k+1 = x_k+u(x_k)
        err = hypot(ux, ux);
        printf("%4d, % .15e, % .15e, % .15e\n", n, *x, *y, err);
    } while (err >= eps);
    return 0;
}
//  ------------    Main    ------------
int main(void)
{
    double x=1, y=0.5, eps = 1e-10;
    char s[128];

    fprintf(stderr,"# x0 y0 = "); fgets(s,128,stdin);
    sscanf(s,"%lf%lf", &x, &y);
    newton2(&x, &y, eps);
    printf("\n# (x, y) = (% .15e, % .15e)\n", x, y);
    return 0;
}
/*  ---   実行例   ---
Z:\proen1\kadai4>newton2_1
# x0 y0 = 1 0.5
# n, x, y, err
   0,  1.000000000000000e+000,  5.000000000000000e-001
   1,  8.500000000000000e-001,  5.500000000000000e-001,  2.121320343559643e-001
   2,  8.266083124196609e-001,  5.634235171696150e-001,  3.308084182290978e-002
   3,  8.260316864653980e-001,  5.636240108887777e-001,  8.154722449349689e-004
   4,  8.260313576542963e-001,  5.636241621612145e-001,  4.650091195410848e-007
   5,  8.260313576541869e-001,  5.636241621612584e-001,  1.546063024937228e-013

# (x, y) = ( 8.260313576541869e-01,  5.636241621612584e-01)
*/