#include <stdio.h>
#include <math.h>
/* Guss-Jordan */
#define N 4
int main() {
int i, j, k, x, y;
double r, max;
double a[N][N + 1] = {
{ 6.0, 2.0, 1.0, 0.0, 26.0 },
{ 3.0, -8.0, -1.0, 3.0, 52.0 },
{ 4.0, 2.0, -9.0, -5.0, 14.0 },
{ 2.0, -9.0, 4.0, -1.0, 78.0 }};
/* forward */
for (i = 0; i < N; i++) {
/* pivotting */
max = fabs(a[i][i]);
k = i;
for (j = i + 1; j < N; j++)
if (max < fabs(a[j][i])) {
max = fabs(a[j][i]);
k = j;
}
if (k != i)
for (j = 0; j <= N; j++) {
r = a[k][j];
a[k][j] = a[i][j];
a[i][j] = r;
}
/* pivotting finished */
r = a[i][i];
for (j = i ; j <= N; j++)
a[i][j] /= r;
for (j = i + 1; j < N; j++) {
r = a[j][i];
for (k = 0; k <= N; k++) {
a[j][k] /= r;
a[j][k] -= a[i][k];
}
}
}
/* backward */
for (i = N - 1; i >= 0; i--) {
for (j = i - 1; j >= 0; j--) {
a[j][N] -= a[i][N] * a[j][i];
a[j][i] = 0.0;
}
}
for (y = 0; y < N; y++) {
for (x = 0; x <= N; x++)
printf("%+.2f ", a[y][x]);
putchar('\n');
}
putchar('\n');
return 0;
}
/* end */