#include <stdio.h>
#include <math.h>
#include <iostream>
#define N    4    // size of matrix
#define N1   N+1
float matrix[N][N1]=
    {{16,-8,-4,0,-8},
  {-8,13,-4,-3,7},
  {-4,-4,9,0,6},
  {0,-3,0,3,-3}
    };
void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%f\n",matrix[i][N]);
  }
}
int main(void)
{
  // Variables declaration
  register char i,j,k;
  float s[N][N],x[N],y[N],d[N],sum;
  // Printing given matrix
  ShowMatrix();
  // Building matrixes S and D
  for (j=0;j<N;j++)
    for (i=0;i<=j;i++)
    {
      sum=matrix[i][j];
      for (k=0;k<i;k++)
        sum-=s[k][i]*d[k]*s[k][j];
      if (i==j)
      {
        d[i]=(sum>0?1:-1);
        s[i][j]=sqrt(fabs(sum));
      }
      else s[i][j]=sum/(d[i]*s[i][i]);
    }
  // Solving equation Gy=b (G: g[i][j]=s[j][i]*d[j])
  for (i=0;i<N;i++)
  {
    y[i]=matrix[i][N];
    for (j=0;j<i;j++) y[i]-=s[j][i]*d[j]*y[j];
    y[i]/=s[i][i]*d[i];
  }
  // Solving equation Sx=y
  for (i=N-1;i>=0;i--)
  {
    x[i]=y[i];
    for (j=i+1;j<N;j++) x[i]-=s[i][j]*x[j];
    x[i]/=s[i][i];
  }
  // Printing solution
  printf("\nSolution:\n");
  for (i=0;i<N;i++)
  printf("x%d=%f\n",i+1,x[i]);
  system("PAUSE");
}