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

/* x^0[1..2]: 3.7  2.8  -->   3.756834008  2.779849593 */
/* x^0[1..2]: 1.4 -1.5  -->   1.373478353 -1.524964836 */
/* x^0[1..2]: 1.0  1.0  -->    -Infinity    -Infinity  */
/* x^0[1..2]: 0.0  0.0  -->   Matrix singular  */

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

  ret[0] = x[0]+3.*log(x[0])-x[1]*x[1];
  ret[1] = 2*x[0]*x[0]-x[0]*x[1]-5.*x[0]+1.;

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

  return ret;
}

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