/*
    newton2.c: 2変数のニュートン法
    x^2 + y^2 = 1, y = x^3 -> (x, y) = (0.826, 0.5636), (-0.826, -0.5636)
 */
#include <stdio.h> // printf, fprintf, gets, sscanf
#include <math.h>  // hypot(hypotenuse): 直角三角形の斜辺を求める関数

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

double J11(double x, double y)
{
    return 2 * x;     //  df_1(x,y)/dx = 2x
}
double J12(double x, double y)
{
    return 2 * y;     //  df_1(x,y)/dy = 2y
}
double J21(double x, double y)
{
    return 3 * x * x; //  df_2(x,y)/dx = 3x^2
}
double J22(double x, double y)
{
    return -1;        //  df_2(x,y)/dy = -1
}

int ucalc(double x, double y, double *ux, double *uy)
/*
    u = -J^{-1}(x)*f(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 eps)
{
    int n=0;
    double x0, y0, ux, uy, err;

    printf("# n, x, y, err\n");
    printf("%4d, % .15e, % .15e\n", n, *x, *y);
    do {
        n++;
        x0 = *x; y0 = *y;
        ucalc(*x, *y, &ux, &uy);
        *x += ux; *y += uy;        // x_{k+1} = x_k + u(x_k)
        err = hypot(*x - x0, *y - y0);
        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 = 1.0e-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
# 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.104876923847481e-004
   4,  8.260313576542963e-001,  5.636241621612145e-001,  3.619393468469258e-007
   5,  8.260313576541869e-001,  5.636241621612584e-001,  1.178637046206069e-013

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