#include <stdio.h>
#include <stdlib.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 = abs(a[i][i]);
    k = i;
    for (j = i + 1; j < N; j++)
      if (max < abs(a[j][i])) {
        max = abs(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 */

+1.00 +0.00 +0.00 +0.00 +6.00 
+0.00 +1.00 +0.00 +0.00 -6.00 
+0.00 +0.00 +1.00 +0.00 +2.00 
-0.00 +0.00 +0.00 +1.00 -4.00