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

#define N_MAX 100


void Triangle(int, double [][N_MAX], double []);
void Solve(int, double [][N_MAX], double [],double []);

main(){
	double A[N_MAX][N_MAX], x[N_MAX], b[N_MAX];
	int n,i,j;

/* Inserimento dei dati (elementi della matrice A, elementi del vettore b) */
	
	do{
		printf("Inserisci il numero di equazioni del sistema \n");
		scanf("%d",&n);
		
		if(n>N_MAX)printf("Il numero di equazioni e' troppo grande!\n");
		
	}while (n>N_MAX);
		
	printf("\nInserisci i coefficienti del sistema (elementi della matrice A)\n");
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			printf("A[%d][%d] = ", i,j);
			scanf("%lf",&A[i][j]);
		}
	}
	
	printf("\nInserisci i termini noti (elementi del vettore b)\n");
	for(i=0;i<n;i++){
	 	printf("b[%d] = ", i);
		scanf("%lf",&b[i]);
	}
	
/* Triangolarizzazione della matrice A mediante un'apposita funzione */
	
	Triangle(n,A,b);
	
/* Soluzione del sistema con la matrice in forma triangolare */
	
	Solve(n,A,b,x);
	
	
   printf("\nLa soluzione ottenuta e'\n");
	for(i=0;i<n;i++){
		printf("x[%d] = %lf\n",i,x[i]);
	}
	
}


void Triangle(int n, double a[N_MAX][N_MAX],double b[N_MAX]){
	int i,j,k;
	double c;
	
	for(i=0;i<n;i++){
/* Divido l'equazione i per l'elemento diagonale a[i][i] della matrice dei coefficienti*/		
		c=a[i][i];
		for(j=i;j<n;j++){
/* Qui sto dividendo gli elementi della riga i (nota bene: sono dentro il doppio ciclo)*/
			a[i][j]=a[i][j]/c;
		}
/* Qui sto dividendo l'elemento i del vettore dei termini noti (nota bene: sono nel ciclo 
 solo su i) */		
		*(b+i)=*(b+i)/c;

/* Sottraggo la riga che ho diviso per a[i][i] dalle altre (nota bene: sono ancora
 nel ciclo su i) */
		
		for(k=i+1;k<n;k++){
			c=a[k][i];
			for(j=i;j<n;j++){
				a[k][j]=a[k][j]-c*a[i][j];
			}
			*(b+k)=*(b+k)-c*(*(b+i));
		}
	}
}
	
	
	
	
	
	
	
	
	
	
void Solve(int n, double a[N_MAX][N_MAX], double b[N_MAX], double x[N_MAX]){
	int i,k;
	double s;
	
	x[n]=b[n]/a[n][n];
	
	for(k=n-1;k>=0;k--){
		s=0.;	
		for(i=k+1;i<n;i++){
			s+=a[k][i]*x[i];
		}
		x[k]=b[k]-s;
	}
	
}