#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "proto.h" 

/* x^0[1..4]:  1. 1. 1. 1. -->  1.0 2.0 3.0 4.0 */
 
double *F(double* x, int n) {
  double *ret = vector(n);
  if ( ret==NULL ) 
    error_exit("Out of memory",__FILE__,__LINE__);

  ret[0] = x[0]*x[1]+x[2]*x[3]-x[0]-3.*x[2]-2.*x[3]+4.;
  ret[1] = x[0]*x[1]-2.*x[0]-x[1]+x[2]+x[3]-5.;
  ret[2] = x[0]*x[0]+x[1]+x[2]+x[3]*x[3]-22.;
  ret[3] = 2*x[0]-x[1]*x[1]*x[1]-x[2]*x[2]+x[3]+11;

  return ret;
}

double **DF(double* x, int n) {
  double **ret = matrix(n,n);
  if ( ret==NULL ) 
    error_exit("Out of memory",__FILE__,__LINE__);

  ret[0][0] = x[1]-1.; ret[0][1] = x[0];    ret[0][2] = x[3]-3.;  ret[0][3] = x[2]-2.;
  ret[1][0] = x[1]-2.; ret[1][1] = x[0]-1.; ret[1][2] = 1.;       ret[1][3] = 1.;
  ret[2][0] = 2.*x[0]; ret[2][1] = 1.;      ret[2][2] = 1.;       ret[2][3] = 2.*x[3];
  ret[3][0] = 2.; ret[3][1] =-2.*x[1]*x[1]; ret[3][2] = -2.*x[2]; ret[3][3] = 1.;

  return ret;
}

int main() {
  int n=4,i;
  double *f, **df, *x, *dx;
  double errx, tol=1.e-12;

  x = vector(n);
  if ( x==NULL ) 
    error_exit("Out of memory",__FILE__,__LINE__);
  printf("x^0[1..%d]: ",n);
  for (i=0; i<n; i++) 
    scanf("%lf",&x[i]); 
  
  do {
    f = F(x,n);
    df = DF(x,n);
    scale(f,-1.0,n);

    dx = gauss(df,f,n);
    xpay(x,1.0,dx,n);
    errx=sqrt(xyT(dx,dx,n));

    for (i=0; i<n; i++) 
      printf("%12.9f%c",x[i],i<n-1?' ':'\n');

    free_vector(f);
    free_matrix(df,n,n);
  } while (errx > tol);
}