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

/* x^0[1..3]:  0.5  0.5  0.5 -->  0.785196933  0.496611393  0.369922831 */
/* x^0[1..3]: -0.5  0.5  0.5 --> -0.785196933  0.496611393  0.369922831 */
/* x^0[1..3]: -1.0 -1.0 -1.0 --> Divergence  */

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[0]+x[1]*x[1]+x[2]*x[2]-1;
  ret[1] = 2.*x[0]*x[0]+x[1]*x[1]-4.*x[2];
  ret[2] = 3.*x[0]*x[0]-4.*x[1]+x[2]*x[2];

  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] = 2.*x[0]; ret[0][1] = 2.*x[1]; ret[0][2] = 2.*x[2]; 
  ret[1][0] = 4.*x[0]; ret[1][1] = 2.*x[1]; ret[1][2] = -4.; 
  ret[2][0] = 6.*x[0]; ret[2][1] = -4.;     ret[2][2] = 2.*x[2]; 

  return ret;
}

int main() {
  int n=3,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);
}