#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define N 3

int main() {
  int A[N][N] = {0}, B[N][2*N] = {0}, b[N] = {0};
  int x[N] = {0};
  /* genera la soluzione, la matrice caratteristica e la colonna dei termini
     noti */
  int i, j;
  for (i = 0; i < N; i++) {
    x[i] = lrand48() % 11;
    printf("x_%d = %d\n", i, x[i]);
    for (j = 0; j < N; j++) {
      A[i][j] = -9 + 2*(lrand48() % 10);
      B[i][j]     = A[i][j];
      B[i][j + N] = A[i][j];
    }
  }
  for (i = 0; i < N; i++) {
    for (j = 0; j < N; j++) {
      b[i] += A[i][j]*x[j]; 
    }
  }
  /* stampa il sistema */
  for (i = 0; i < N; i++) {
    for (j = 0; j < N; j++) {
      if (A[i][j] != 0) {
	if (abs(A[i][j]) != 1) {
	  printf("%+dx_%d ", A[i][j], j);
	} else {
	  if (A[i][j] > 0) {
	    printf("+");
	  } else {
	    printf("-");
	  }
	  printf(" x_%d ", j);
	}
      } else {
	printf("      ");
      }
    }
    printf(" = %+d\n", b[i]);
  }
  /* applica la regola di Sarrus per il calcolo del determinante di A */
  int detA = 0;
  for (i = 0; i < N; i++) {
    int prodp = 1;
    int prodm = 1;
    for (j = 0; j < N; j++) {
      prodp *= B[j][j + i];
      prodm *= B[j][N + 1 - j - i];
    }
    detA += prodp - prodm;
  }
  printf("\ndet(A) = %d\n", detA);
  /* applica la regola di Sarrus alla matrice con la colonna j sostituita */
  int detAj[N] = {0};
  int C[N][2*N];
  int k;
  for (k = 0; k < N; k++) {
    for (i = 0; i < N; i++) {
      for (j = 0; j < 2*N; j++) {
	C[i][j] = B[i][j];
	if ((k == j) || (k == (j - N))) {
	  C[i][j] = b[i];
	}
      }
    }
    for (i = 0; i < N; i++) {
      int prodp = 1;
      int prodm = 1;
      for (j = 0; j < N; j++) {
	prodp *= C[j][j + i];
	prodm *= C[j][N + 1 - j - i];
      }
      detAj[k] += prodp - prodm;
    }
    printf("\ndet(A_%d) = %d", k, detAj[k]);
  }
  printf("\n");
  for (i = 0; i < N; i++) {
    float y = (float)detAj[i]/detA;
    printf("\nx_%d = %f (%d) diff=%f err=%f", i, y, x[i], fabs(x[i]-y), fabs(x[i]-y)/x[i]);
  }
  printf("\n");
}