/*
    newton2_3.c: 関数ポインタ追加バージョン
    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 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);  // 配列f[]に f(x) の各成分を格納
    Jcalc(x, J); // 配列J[][]に J(x) の各成分を格納
    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]; // x_k+1 = x_k+u(x_k)
        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, func, Jcalc, eps);
    printf("\n# (x, y) = (% .15e, % .15e)\n", x[1], x[2]);
    return 0;
}
/*  ---   実行例   ---
Z:\proen1\kadai4>newton2_3
# x0 y0 = 1 0.5
# n, x, y, err
   0,  1.000000000000000e+000,  5.000000000000000e-001
   1,  8.500000000000000e-001,  5.500000000000000e-001,  1.581138830084190e-001
   2,  8.266083124196609e-001,  5.634235171696150e-001,  2.696964703252782e-002
   3,  8.260316864653980e-001,  5.636240108887777e-001,  6.104876923847893e-004
   4,  8.260313576542963e-001,  5.636241621612145e-001,  3.619393467961632e-007
   5,  8.260313576541869e-001,  5.636241621612584e-001,  1.178466214256035e-013

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