#include <stdio.h>
#include <math.h>
/* gauss.c : 1999 5 8 by Oguni */
#define N 3
main(){
  double A[N][N]={1.,4.,3.,2.,5.,4.,1.,-3.,-2.}, Aa[N][N];
  double b[N]={1.,4.,5.}, x[N], bb[N], e[N];
  double tol=0.000001, s;
  int n=3, i, j;
  for (i=0; i<n; i++){
    for (j=0; j<n; j++){
      Aa[i][j]=A[i][j]; }
    bb[i]=b[i];
  }
  gauss(A, b, x, tol, n);
  for (i=0; i<n; i++){
    s=0.0;
    for (j=0; j<n; j++){
      s+=Aa[i][j]*x[j];}
    e[i]=s-bb[i];
    printf("x(%d)=%f error=%f\n",i, x[i], e[i]); }
}
/* function */
gauss(A, b, x, tol, n)
double A[N][N], b[N], x[N], tol;
int n;
{
  int i, j, k;
  double akk, aik, s;
  for (k=0; k<n-1; k++){
    if (fabs(A[k][k])<tol){
      printf("The pivot is too small at step # %d", k); return; }
    akk=1.0/A[k][k];
    for (i=k+1; i<n; i++){
      aik=-A[i][k]*akk;
      for (j=k+1; j<n; j++){
        A[i][j]+=aik*A[k][j]; }
      b[i]+=aik*b[k];
    }
    for (j=k+1; j<n; j++){
      A[k][j]*=akk; }
    b[k]*=akk;
  }
  x[n-1]=b[n-1]/A[n-1][n-1];
  for (k=n-2; k>=0; k--){
    s=0.0;
    for (j=k+1; j<n; j++){ s+=A[k][j]*x[j]; }
    x[k]=b[k]-s;
  }
}