/* 連立非線型方程式のNewton 法*/
/* sim newton2.c */
#include <stdio.h>
#include <math.h>
#define EPS 1.0e-12
#define NUMAX 10
#define N 2
void gauss(double (*a)[N],double *b);
double f(int i, double x,double y);
double df(int i,int j,double x,double y);

int main(void)
{
  int i,j,nu;
  double jacobi[N][N];
  double fx[N];
  double x,y;
  double err,err1,err2;

  nu = 0;
  x = -1.0;//1.0;//4.0;
  y = -1.0;//1.0;//3.0;
  err1 = fabs(f(0,x,y));
  err2 = fabs(f(1,x,y));
  if (err1 > err2)
    err = err1;
  else
    err = err2;
  printf("%3d %20.15f %20.15f %10.3e\n",nu,x,y,err);
  while(err > EPS && nu < NUMAX){
    for(i = 0;i < N;i++){
      for(j = 0;j < N;j++){
	jacobi[i][j] = df(i,j,x,y);
      }
      fx[i] = -f(i,x,y);
    }
    gauss(jacobi,fx);
    x = fx[0] + x;
    y = fx[1] + y;
    nu++;
    err1 = fabs(f(0,x,y));
    err2 = fabs(f(1,x,y));
    if (err1 > err2)
      err = err1;
    else
      err = err2;
    printf("%3d %20.15f %20.15f %10.3e\n",nu,x,y,err);
  }
  return 0;
}

void gauss(double (*a)[N],double *b)
{
  double mul,w;
  int i,j,k,p;
  for(k = 0; k < N - 1; k++){
    p = k;
    for(i = k + 1; i < N; i++){
      if (fabs(a[p][k]) < fabs(a[i][k])) p = i;
    }
    if (p > k){
      for(j = k; j < N; j++){
	w = a[p][j];
	a[p][j] = a[k][j];
	a[k][j] = w;
      }
      w = b[p];
      b[p] = b[k];
      b[k] = w;
    }
    for(i = k + 1; i < N; i++){
      mul = a[i][k]/a[k][k];
      for(j = k; j < N; j++)
	a[i][j] -= mul * a[k][j];
      b[i] -= mul * b[k];
    }
  }
  b[N-1] /= a[N-1][N-1];
  for(i = N - 2; i>= 0; i--){
    w = 0.0;
    for(k = i + 1;k < N;k++)
      w += a[i][k] * b[k];
    b[i] = (b[i] - w) / a[i][i];
  }
}

double f(int i,double x,double y)
{
  double z;
  switch(i){
  case 0:
    //    z = 2 * x - y * y + log(x);
    z = x * x - 2 * x * y + 3 * y * y -1;
    break;
  case 1:
    //    z = x * x - x * y - x + 1.0;
    z = x * exp(x) + y * y - 3 * x * y -1;
  }
  return z;
}

double df(int i,int j,double x,double y)
{
  double z;
  int k;
  k = N * i + j;
  switch(k){
  case 0:
    //    z = 2 + 1.0 /x;
    z = 2 * x - 2 * y;
    break;
  case 1:
    //    z = - 2.0 * y;
    z = -2 * x + 6 * y;
    break;
  case 2:
    //    z = 2.0 * x - y - 1.0;
    z = exp(x) + x * exp(x) -3 * y;
    break;
  case 3:
    //    z = - x;
    z = 2 * y - 3 * x;
  }
  return z;
}