#include <stdio.h>
#include <math.h>

#define N_STOP 100
#define TOL 1.0e-16


double g_xprime(double x, double y);
double g_yprime(double x, double y);
double g_xxprime(double x, double y);
double g_xyprime(double x, double y);
double g_yyprime(double x, double y);

int main ()
{
    int i=0;
    double a,b,c,det,x,y,x_next,y_next;
    x=0.1;
    y=1.3;
   
    while((i < N_STOP) && ((fabs(x-0.) > TOL) || (fabs(y-1.0) > TOL)) )
    {
                         //compute the second derivatives of g at (x,y)
       a=g_xxprime(x,y);
       b=g_yyprime(x,y);
       c=g_xyprime(x,y);

       det=a*b-c*c;	 //compute the determinant of Hessian
       
                         //next iteration
       x_next=x-(b*g_xprime(x,y)-c*g_yprime(x,y))/det;
       y_next=y-(-c*g_xprime(x,y)+a*g_yprime(x,y))/det;
       
       i++;
       x=x_next;
       y=y_next;
       printf("x_%d =%.16f   y_%d=%.16f \n",i,x,i,y);
    } 

    return 0;
}


double g_xprime(double x,double y)
{
     return -2.*x*exp(-x*x-y*y)*(x*x+2.*y*y-1.);
}

double g_yprime(double x, double y)
{
     return -2.*y*exp(-x*x-y*y)*(x*x+2.*y*y-2.);
}

double g_xxprime(double x, double y)
{
     return 2.*exp(-x*x-y*y)*(-5.*x*x-2.*y*y+2.*x*x*x*x+4.*x*x*y*y+1.);
}

double g_xyprime(double x, double y)
{
     return 4.*x*y*exp(-x*x-y*y)*(x*x+2.*y*y-3.);
}

double g_yyprime(double x, double y)
{
     return 2.*exp(-x*x-y*y)*(-x*x-10.*y*y+2.*x*x*y*y+4.*y*y*y*y+2.);
}