/* Gaussian Elimination */
#include <math.h>
#include <stdio.h>

main()
{
double matrix[16][16];
double b[16],x[16];
int n,i,count,j;
double ratio,temp;

/* Read in matrix and b */
scanf("%d",&n);
printf("The matrix is %d by %d\n",n,n);
printf("Matrix:\n");
for(i=0;i<n;i++){
        for(j=0;j<n;j++){
                /*scanf("%lf ",&matrix[i][j]);*/
                /* for Hilbert matrix */
                matrix[i][j] = 1. / (i+j+1);
        }
        scanf("%lf ",&b[i]);
        scanf("\n");
        b[i] = 1.;
}
for(i=0;i<n;i++){
        for(j=0;j<n;j++){
                printf("%lf ",matrix[i][j]);
        }
        printf("\n");
}
printf("Right Hand Side:\n");
for(i=0;i<n;i++){
        printf("%lf ",b[i]);
}
printf("\n");
/* Gaussian elimination */
for(i=0;i<(n-1);i++){
        for(j=(i+1);j<n;j++){
                ratio = matrix[j][i] / matrix[i][i];
                for(count=i;count<n;count++) {
                        matrix[j][count] -= (ratio * matrix[i][count]);
                }
                b[j] -= (ratio * b[i]);
        }
}
for (i=0;i<=n-1;i++){
        for(j=0;j<n;j++){
                printf("%lf ",matrix[i][j]);
        }
        printf("\n");
}
/* Back substitution */
x[n-1] = b[n-1] / matrix[n-1][n-1];
for(i=(n-2);i>=0;i--){
        temp = b[i];
        for(j=(i+1);j<n;j++){
                temp -= (matrix[i][j] * x[j]);
        }
        x[i] = temp / matrix[i][i];
}
printf("Answer:\n");
for(i=0;i<n;i++){
        printf("x%d = %lf\n",i,x[i]);
}
}