/*求函數y=x1^2+3*x2^2-4x1的極小值*/

#include"stdio.h"
#include"math.h"
#define N 50

/*a[N]表示梯度，b[N][N]表示Hesse矩陣*/
void function(float x1,float x2,float a 
[N],float b[N][N])
{
 a[0]=2*x1-4;
 a[1]=6*x2;
 b[0][0]=2;
 b[0][1]=0;
 b[1][0]=0;
 b[1][1]=6;
}

/*求矩陣的逆*/
void Reverse(float a[N][N],float c[N] 
[N],int n)
{int i,j,k,t,m,q,r,h;
 float b[N][2*N];
 float p=0.0;
 for(i=0;i<n;i++)
  for(j=0;j<n;j++)
   b[i][j]=a[i][j];
 for(i=0;i<n;i++)
  for(j=n;j<2*n;j++)
   {if(i==j-n)b[i][j]=1.0;
   else b[i][j]=0.0;}
 for(t=0;t<n;t++)
   {
    for(k=t+1;k<n;k++)
    for(j=t+1;j<2*n;j++)
     {if(b[t][t]==0)
     	{for(k=t;k<n;k++)
     	 if(b[k][t]!=0){r=k;k=t+1;break;}
     	 for(m=0;m<2*n;m++)
     	{p=b[r][m];
     	 b[r][m]=b[t][m];
     	 b[t][m]=p;}
         for(i=0;i<2*n;i++)
          if(i!=t)
     	  b[t][i]=b[t][i]/b[t][t];
     	  b[t][t]=1;
     	b[k][j]=b[k][j]-b[t][j]*b[k][t];
     	}
      else { 	
        for(i=0;i<2*n;i++)
    	  if(i!=t)
     	  b[t][i]=b[t][i]/b[t][t];
          b[t][t]=1;
     	b[k][j]=b[k][j]-b[t][j]*b[k][t];
     	}
     }
     }
   for(i=0;i<n;i++)
   for(j=i+1;j<n;j++)
   b[j][i]=0;
   for(j=n;j<2*n;j++)
   b[n-1][j]=b[n-1][j]/b[n-1][n-1];
   b[n-1][n-1]=1;
   for(j=n-1;j>=1;j--)
   for(i=j-1;i>=0;i--)
   for(k=n;k<2*n;k++)
   b[i][k]=b[i][k]-b[j][k]*b[i][j];
   for(j=0;j<n;j++)
   for(i=j+1;i<n;i++)
   b[j][i]=0;
   for(i=0;i<n;i++)
   for(j=n;j<2*n;j++)
    c[i][j-n]=b[i][j];
 }
 

/*牛頓法*/
void Newton(float t)/*t表示精度*/
{float t1[2],t2[2],sum=0; 
 float a[N], b[N][N],c[N][N];
 int i,j;
 printf("Please enter the first point:\n");
 for(i=0;i<2;i++)
 scanf("%f",&t1[i]);/*輸入初始點*/
 while(1)
  { function(t1[0],t1[1],a,b);
    Reverse(b,c,2);
  	for(i=0;i<2;i++)
   { sum=0;
    for(j=0;j<2;j++)
     sum+=c[i][j]*a[j];
    t2[i]=t1[i]-sum;
    }
     sum=0;
     sum=sqrt(pow(fabs(t2[0]-t1[0]),2)+pow(fabs(t2[1]-t1[1]),2));
    if(sum<t)
     {
      printf("the minimum point is:\n");  
      for(i=0;i<2;i++)
      printf("%f ",t2[i]);
      break;
     }
    else 
     {for(i=0;i<2;i++)
      t1[i]=t2[i];
     }
  }
}

 void main()
{float t;/*t表示精度*/
 printf("Please enter the precision:\n");
 scanf("%f",&t);
 Newton(t);
}
/*注：初始值(x1,x2),x2不能為0*/