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


typedef struct { double val[3]; } vect;


typedef double (*Func)(vect);

#define EPSILON 0.00001
#define N 3                    //未知数の数

int printA(double a[N][N+1]);
void gauss(double a[N][N+1], double x[N]);
int gauss_jordan(int n, double a[N][N+1], double b[N]);


int printV(vect X){ for(int i=0;i<3;i++) printf("%10.3f",X.val[i]);	  printf("\n"); }


double f1(vect X){ 
double x=X.val[0],y=X.val[1],z=X.val[2];
return -y*y * z*z - y*y + 24*y*z - z*z -13; }

double f2(vect X){ 
double x=X.val[0],y=X.val[1],z=X.val[2];
return -x*x * z*z - x*x + 24*x*z - z*z -13; }

double f3(vect X){ 
double x=X.val[0],y=X.val[1],z=X.val[2];
return -x*x * y*y - x*x + 24*x*y - y*y -13; }


vect d2v(double a,double b,double c) { vect x; x.val[0]=a; x.val[1]=b; x.val[2]=c; return x; }
vect vect_add(vect x, vect y) { vect z; z.val[0]=x.val[0]+y.val[0];  z.val[1]=x.val[1]+y.val[1];  z.val[2]=x.val[2]+y.val[2]; return z; }
vect vect_hiku(vect x, vect y) { vect z; z.val[0]=x.val[0]-y.val[0];  z.val[1]=x.val[1]-y.val[1];  z.val[2]=x.val[2]-y.val[2]; return z; }

double bibun_x(Func f, vect a) { return  ( f(vect_add(a,d2v(EPSILON,0,0))) - f(a) )/EPSILON;  }
double bibun_y(Func f, vect a) { return  ( f(vect_add(a,d2v(0,EPSILON,0))) - f(a) )/EPSILON;  }
double bibun_z(Func f, vect a) { return  ( f(vect_add(a,d2v(0,0,EPSILON))) - f(a) )/EPSILON;  }



 
int main() {
vect V=d2v(1,1,1);
double a[3][4], X[3];

    for (int i = 0; i < 10; i++) {
  a[0][0] = bibun_x(f1,V);  a[0][1] = bibun_y(f1,V);  a[0][2] = bibun_z(f1,V);
  a[1][0] = bibun_x(f2,V);  a[1][1] = bibun_x(f2,V);  a[1][2] = bibun_x(f2,V);
  a[2][0] = bibun_x(f3,V);  a[2][1] = bibun_x(f3,V);  a[2][2] = bibun_x(f3,V);
  a[0][3] = V.val[0];   a[1][3] = V.val[1];   a[2][3] = V.val[2];
  X[0] = V.val[0];   X[1] = V.val[1];   X[2] = V.val[2];
//printV(V);printA(a);gauss(a,X);

 gauss_jordan(3, a, X);

V = vect_hiku(V, d2v(X[0],X[1],X[2]));
printV(V);
}
return 0;
}




//http://s...content-available-to-author-only...e.jp/wiki/wiki.cgi?page=%CF%A2%CE%A9%B0%EC%BC%A1%CA%FD%C4%F8%BC%B0%A4%CE%B2%F2%CB%A1


void gauss(double a[N][N+1], double x[N]){
        int i,j,k,l,pivot;
        double p,q,m,tmp;

	/*縦方向のループ*/
        for(i=0;i<N;i++) {
	  m=0;
	  pivot=i;
	  
	  /*i列中でi行目より下にある要素内で最大を選ぶ*/
	  for(l=i;l<N;l++) {
	    /*i列の中で一番値が大きい行を選ぶ*/
	    if(fabs(a[l][i])>m) {
	      m=fabs(a[l][i]);
	      pivot=l;
	    }
	  }

	  /*pivotがiと違えば、行の入れ替え(これを部分ピボットと呼ぶ)*/
	  if(pivot!=i) {  
	    for(j=0;j<N+1;j++) {
	      tmp=a[i][j];        
	      a[i][j]=a[pivot][j];
	      a[pivot][j]=tmp;
	    }
	  }
        }
	
        for(k=0;k<N;k++) {
	  p=a[k][k];              //対角要素を保存
	  /*対角要素は1になることがわかっているので直接代入*/
	  a[k][k]=1;      
	  
	  /*k行目でk+1列より右にあるすべての要素を対角成分で割る*/
	  for(j=k+1;j<N+1;j++) {
	    a[k][j]/=p;
	  }
	  
	  for(i=k+1;i<N;i++) {
	    q=a[i][k];
	    
	    for(j=k+1;j<N+1;j++) {
	      a[i][j]-=q*a[k][j];
	    }
	    /*0となることがわかっているので直接代入*/
	    a[i][k]=0;
	  }
        }
	
	/*解の計算*/
        for(i=N-1;i>=0;i--) {
	  x[i]=a[i][N];
	  for(j=N-1;j>i;j--) {
	    x[i]-=a[i][j]*x[j];
	  }
        }

}


int printA(double a[N][N+1]){
        int i,j;
        for(i=0;i<N;i++) {
	  for(j=0;j<N+1;j++) {
	    printf("%10.3f",a[i][j]);
	  }
	  printf("\n");
	  
        }
}





//http://w...content-available-to-author-only...0.jp/yamamoto/lecture/2006/5E/Linear_eauations/gaussj_html/node2.html#SECTION00021000000000000000
int gauss_jordan(int n, double a[N][N+1], double b[N]) {
	   int ipv, i, j;
	   double inv_pivot, temp;
	   double big;
	   int pivot_row, row[30];

	   for(ipv=1 ; ipv <= n ; ipv++){

	      /* ---- 最大値探索 ---------------------------- */
	      big=0.0;
	      for(i=ipv ; i<=n ; i++){
	         if(fabs(a[i][ipv]) > big){
	            big = fabs(a[i][ipv]);
	            pivot_row = i;
	         }
	      }
	      if(big == 0.0){return 0;}	
	      row[ipv] = pivot_row;

	      /* ---- 行の入れ替え -------------------------- */
	      if(ipv != pivot_row){
	         for(i=1 ; i<=n ; i++){
	            temp = a[ipv][i];
	            a[ipv][i] = a[pivot_row][i];
	            a[pivot_row][i] = temp;
	         }
	         temp = b[ipv];
	         b[ipv] = b[pivot_row];
	         b[pivot_row] = temp;
	      }

	      /* ---- 対角成分=1(ピボット行の処理) ---------- */
	      inv_pivot = 1.0/a[ipv][ipv];
	      a[ipv][ipv]=1.0;
	      for(j=1 ; j <= n ; j++){
	         a[ipv][j] *= inv_pivot;
	      }
	      b[ipv] *= inv_pivot;

	      /* ---- ピボット列=0(ピボット行以外の処理) ---- */
	      for(i=1 ; i<=n ; i++){
	         if(i != ipv){
	            temp = a[i][ipv];
		    a[i][ipv]=0.0;
	            for(j=1 ; j<=n ; j++){
	               a[i][j] -= temp*a[ipv][j];
	            }
	            b[i] -= temp*b[ipv];
	         }
	      }

	   }

	   /* ---- 列の入れ替え(逆行列) -------------------------- */
	   for(j=n ; j>=1 ; j--){
	      if(j != row[j]){
	         for(i=1 ; i<=n ; i++){
	            temp = a[i][j];
	            a[i][j]=a[i][row[j]];
	            a[i][row[j]]=temp;
	         }
	      }
	   }

	   return 1;
	}