/* alg_lin_simpl.c : inizialmente, si determina la matrice Q associata
                     all'approssimazione quadratica della Hamiltoniana
		     che descrive il modello ristretto, circolare, piano
		     dei tre corpi, nell'intorno di uno dei due punti
		     di equilibrio equilateri L4 o L5.
		     Successivamente, per tappe successive, si costruisce
		     la matrice simplettica U che pone Q in forma
		     canonica (cioe' quella associata a una Hamiltoniana
		     quadratica che descrive una coppia di oscillatori
		     armonici).
*/

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

/* E' necessario includere il file stdlib.h, al fine di utilizzare la
   function exit che provoca la fine immediata dell'esecuzione del
   programma */
#include <stdlib.h>

/* La prossima riga fa si' che ogni volta che viene riportata la
   scritta NDIM all'interno del programma, questa stessa scritta in
   realta' viene sostituita dal valore riportato accanto (in questo
   caso 4), per effetto dell'azione del preprocessore */
#define NDIM 4


/* Definizione dei prototipi delle funzioni */
void gramschmidt_simpl_duevet(double v1[], double v2[], int n);
void inv_mat_cramer(double A[][NDIM], int n, double A_inv[][NDIM]);
double determinante(double mat[][NDIM], int n);
double met_potenze_simpl(double A[][NDIM], int n, double v1[], double v2[]);
double norma_euclidea(double v[], int n);
void prod_matr(double A[][NDIM], double B[][NDIM], int n, double C[][NDIM]);
void prod_matr_vet(double mat[][NDIM], double vet[], int n, double vetris[]);
double prod_scal(double v1[], double v2[], int n);
int risolvi_eq(double a, double b, double c, double *x1, double *x2);
void togli_riga0_colonnai(double mat[][NDIM], int n,
			  int i, double matridotta[][NDIM]);
void togli_rigal_colonnal(double mat[][NDIM], int n,
			  int l, double matridotta[][NDIM]);
void trasp_matr(double A[][NDIM], int n, double At[][NDIM]);


/* ### */
int main()
{


  printf("\n     ###   Obiettivo 1   ###     \n");
  int i, j, n=NDIM/2, risp;
  double mu, gamma;
  double Q[NDIM][NDIM], J[NDIM][NDIM], JQ[NDIM][NDIM];

  /* Input del valore del parametro mu */
  do {
    printf("Valore del rapporto delle masse mu? (con 0 < mu < ");
    printf("[1-sqrt(69)/9]/2 )\n    ");
    scanf("%lf", &mu);
  } while(mu <= 0. || mu >= (1-sqrt(69.)/9.)/2. );
  /* Scelta del punto di equilibrio Lagrangiano */
  do {
    printf("Studio dell'intorno del punto L4 oppure di quello di L5? ");
    printf("(si risponda 4 o 5)\n    ");
    scanf("%d", &risp);
  } while(risp < 4 || risp >5);
  /* Definizione del parametro gamma */
  gamma = -3. * sqrt(3.) / 4. * (1. - 2. * mu);
  if(risp == 5)
    gamma *= -1.;
  /* Definizione della matrice Q */
  Q[0][0] = 1.;  Q[0][1] = 0.;  Q[0][2] = 0.;  Q[0][3] = 1.;  
  Q[1][0] = 0.;  Q[1][1] = 1.;  Q[1][2] = -1.;  Q[1][3] = 0.;  
  Q[2][0] = 0.;  Q[2][1] = -1.;  Q[2][2] = 0.25;  Q[2][3] = gamma;  
  Q[3][0] = 1.;  Q[3][1] = 0.;  Q[3][2] = gamma;  Q[3][3] = -1.25;  
  /* Definiamo J uguale alla matrice fondamentale della geometria
     simplettica */
  for(i=0; i<2*n; i++)
    for(j=0; j<2*n; j++)
      J[i][j] = 0.;
  for(i=0; i<n; i++)
    J[i][i+n] = -1;
  for(i=0; i<n; i++)
    J[i+n][i] = 1;
  /* Eseguiamo il prodotto matriciale JQ = J * Q, facendo ricorso alla
     apposita function prod_matr */
  prod_matr(J, Q, 2*n, JQ);
  /* Stampa (ordinata) della matrice JQ */
  printf("  Stampa della matrice JQ:\n");
  for(i=0; i<2*n; i++) {
    printf("(");
    for(j=0; j<2*n; j++)
      printf(" %20.16lf", JQ[i][j]);
    printf(" )\n");
  }




  printf("\n\n\n     ###   Obiettivo 2   ###     \n");
  double alpha[5];
  double mat1[NDIM][NDIM], mat2[NDIM][NDIM];
  /*    ### Calcolo dei coefficienti alpha[i], associati
	    al polinomio caratteristico della matrice JQ  ###       */
  /* Calcolo del coefficiente alpha[0], che deve essere posto uguale
     al determinante della matrice JQ */
  alpha[0] = determinante(JQ, NDIM);
  /* Calcolo del coefficiente alpha[1], che viene inizialmente posto
     uguale a 0 */
  alpha[1] = 0;
  for(i=0;i<NDIM;i++) {
    /* Poniamo mat1 uguale alla matrice che si ottiene da JQ togliendo
       la i-esima riga e i-esima colonna. Cio' viene effettuato
       facendo ricorso alla function togli_rigal_colonnal */
    togli_rigal_colonnal(JQ, NDIM, i, mat1);
    /* Ridefiniamo alpha[1] sottraendo al suo valore precedente il
       determinante di mat1 */ 
    alpha[1] -= determinante(mat1, NDIM-1);
  }
  /* Calcolo del coefficiente alpha[2], che viene inizialmente posto
     uguale a 0 */
  alpha[2] = 0;
  for(i=0;i<NDIM;i++) {
    /* Poniamo mat1 uguale alla matrice che si ottiene da JQ togliendo
       la i-esima riga e i-esima colonna. Cio' viene effettuato
       facendo ricorso alla function togli_rigal_colonnal */
    togli_rigal_colonnal(JQ, NDIM, i, mat1);
    for(j=i;j<NDIM-1;j++) {
      /* Poniamo mat2 uguale alla matrice che si ottiene da mat1 togliendo
	 la j-esima riga e j-esima colonna. Cio' viene effettuato
	 facendo ricorso alla function togli_rigal_colonnal */
      togli_rigal_colonnal(mat1, NDIM-1, j, mat2);
      /* Ridefiniamo alpha[2] sommando al suo valore precedente il
	 determinante di mat2 */ 
      alpha[2] += determinante(mat2, NDIM-2);
    }
  }
  /* Calcolo del coefficiente alpha[3], che deve essere posto uguale a
     meno la traccia della matrice A */
  alpha[3] = 0;
  for(i=0;i<NDIM;i++)
    alpha[3] -= JQ[i][i];
  /* Poniamo il coefficiente alpha[4] = 1 */
  alpha[4] = 1;

  /* Stampa dei coefficienti alpha[i], associati al polinomio
     caratteristico della matrice JQ */
  printf("  Stampa dei coefficienti alpha[i], ciascuno dei quali e' associato\n");
  printf("  al termine di i-esimo grado del polinomi caratteristico della matrice JQ:\n");
  for(i=NDIM;i>=0;i--) {
    printf("    alpha[%d] = %le\n", i, alpha[i]);
  }




  printf("\n\n\n     ###   Obiettivo 3   ###     \n");
  int flag;
  double a, b, c, x1, x2;
  double omega[NDIM/2];
  /* Definizione dei coefficienti a, b, c dell'equazione di secondo
     grado associata al polinomio caratteristico (che e' biquadratico)
     della matrice JQ */
  a = alpha[4];
  b = alpha[2];
  c = alpha[0];
  /* Risolviamo l'equazione di secondo grado, assegnando alle
     variabili x1 e x2 le soluzioni e a flag il valore 1, 0, -1
     a seconda del segno di b^2 - 4*a*c */
  flag = risolvi_eq(a, b, c, &x1, &x2);
  if(flag == -1) {
    printf("  Non e' possibile che esistano le seguenti due soluzioni");
    /* Stampiamo le soluzioni nella forma x1 +- i*x2 */
    printf(" COMPLESSE\n");
    printf(" dell'equazione di secondo grado:    x1,x2 = %lf +- i*%lf\n",
	   x1, x2);
    /* Arresto forzato dell'esecuzione del programma */
    exit(1);
  }
  if( x1 > 0 || x2 > 0 ) {
    printf("  Non e' possibile che una delle seguenti soluzioni dell'equazione\n");
    /* Stampiamo le soluzioni nella forma x1 +- i*x2 */
    printf("  di secondo grado sia positiva:\n    x1 = %lf    x2 = %lf\n",
	   x1, x2);
    /* Arresto forzato dell'esecuzione del programma */
    exit(1);    
  }
  printf("    Gli autovalori della matrice JQ sono:\n");
  printf("      +-i * %lf    e    -+i * %lf\n", sqrt(-x1), sqrt(-x2));
  omega[0] = sqrt(-x1);
  omega[1] = -sqrt(-x2);


  printf("\n\n\n     ###   Obiettivo 4   ###     \n");
  double JQJQ[NDIM][NDIM], v1[NDIM], v2[NDIM], vtmp[NDIM], lambda_barrato;

  /* Eseguiamo il prodotto matriciale JQJQ = JQ * JQ, facendo ricorso
     alla apposita function prod_matr */
  prod_matr(JQ, JQ, 2*n, JQJQ);
  /* Calcoliamo l'autovalore lambda_barrato che e' quello massimo (in
     valore assoluto) della matrice JQJQ e 2 autovettori v1 e v2 che
     appartengono all'autospazio ad esso associato. v1 e v2 sono tali
     che la norma euclidea di v1 e' uguale a 1 e il prodotto scalare
     tra v1 e J * v2 e' uguale a -1. Tutto questo calcolo viene
     effettuato facendo ricorso alla apposita function
     met_potenze_simpl */
  lambda_barrato = met_potenze_simpl(JQJQ, 2*n, v1, v2);
  printf("  Si ricordi che l'autovalore lambda_barrato che e' massimo (in\n");
  printf("  valore assoluto) per la matrice JQJQ deve essere uguale al\n");
  printf("  quadrato di omega[0] cambiato di segno; noi abbiamo ottenuto che:\n");
  printf("    lambda_barrato = %lf\n", lambda_barrato);
  /* Definisco il vettore vtmp in modo tale che vtmp = JQJQ * v1 */
  prod_matr_vet(JQJQ, v1, 2*n, vtmp);
  /* Sommiamo a vtmp il vettore che si ottiene moltiplicando scalarmente
     il quadrato di omega[0] per l'autovettore v1 */
  for(i=0;i<2*n;i++)
    vtmp[i] += omega[0] * omega[0] * v1[i];
  printf("  Verifichiamo che v1 e' effettivamente un autovettore della matrice\n");
  printf("  JQJQ per l'autovalore -omega[0]*omega[0] a meno del seguente errore:\n");
  printf("    %le\n", norma_euclidea(vtmp, 2*n) );
  /* Definisco il vettore vtmp in modo tale che vtmp = JQJQ * v2 */
  prod_matr_vet(JQJQ, v2, 2*n, vtmp);
  /* Sommiamo a vtmp il vettore che si ottiene moltiplicando scalarmente
     il quadrato di omega[0] per l'autovettore v2 */
  for(i=0;i<2*n;i++)
    vtmp[i] += omega[0] * omega[0] * v2[i];
  printf("  Verifichiamo che v2 e' effettivamente un autovettore della matrice\n");
  printf("  JQJQ per l'autovalore -omega[0]*omega[0] a meno del seguente errore:\n");
  printf("    %le\n", norma_euclidea(vtmp, 2*n) );
  /* Definisco il vettore vtmp in modo tale che vtmp = J v2 */
  prod_matr_vet(J, v2, 2*n, vtmp);
  printf("  Verifichiamo che il prodotto scalare tra l'autovettore v1 e \n");
  printf("  il vettore J * v2 e' uguale a -1 a meno del seguente errore:\n");
  printf("    %le\n", fabs( prod_scal(v1, vtmp, 2*n) + 1 ) );




  printf("\n\n\n     ###   Obiettivo 5   ###     \n");
  double  S[NDIM][NDIM], inv_JQJQ[NDIM][NDIM], V[NDIM][NDIM], JV[NDIM][NDIM];
  double VtJV[NDIM][NDIM], VtJV_men_J[NDIM][NDIM], max_VtJV_men_J;

  /* Poniamo la prima e la terza colonna di V uguali, rispettivamente
     a v1 e v2 */
  for(i=0;i<2*n;i++) {
    V[i][0] = v1[i];
    V[i][2] = v2[i];
  }
  /* Poniamo inv_JQJQ uguale alla matrice inversa di JQJQ = JQ * JQ,
     facendo ricorso alla function inv_mat_cramer, che e' basata sul
     metodo di Cramer */
  inv_mat_cramer(JQJQ, 2*n, inv_JQJQ);
  /* Calcoliamo l'autovalore lambda_barrato che e' quello massimo (in
     valore assoluto) della matrice inv_JQJQ e 2 autovettori v1 e v2
     che appartengono all'autospazio ad esso associato. v1 e v2 sono
     tali che la norma euclidea di v1 e' uguale a 1 e il prodotto
     scalare tra v1 e J * v2 e' uguale a -1. Tutto questo calcolo
     viene effettuato facendo ricorso alla apposita function
     met_potenze_simpl */
  lambda_barrato = met_potenze_simpl(inv_JQJQ, 2*n, v1, v2);

  printf("  Si ricordi che l'autovalore lambda_barrato che e' massimo (in\n");
  printf("  valore assoluto) per la matrice inv_JQJQ deve essere uguale a\n");
  printf("  -1 fratto il quadrato di omega[1]; noi abbiamo ottenuto che:\n");
  printf("    lambda_barrato = %lf\n", lambda_barrato);
  /* Definisco il vettore vtmp in modo tale che vtmp = JQJQ * v1 */
  prod_matr_vet(JQJQ, v1, 2*n, vtmp);
  /* Sommiamo a vtmp il vettore che si ottiene moltiplicando scalarmente
     il quadrato di omega[1] per l'autovettore v1 */
  for(i=0;i<2*n;i++)
    vtmp[i] += omega[1] * omega[1] * v1[i];
  printf("  Verifichiamo che v1 e' effettivamente un autovettore della matrice\n");
  printf("  JQJQ per l'autovalore -omega[1]*omega[1] a meno del seguente errore:\n");
  printf("    %le\n", norma_euclidea(vtmp, 2*n) );
  /* Definisco il vettore vtmp in modo tale che vtmp = JQJQ * v2 */
  prod_matr_vet(JQJQ, v2, 2*n, vtmp);
  /* Sommiamo a vtmp il vettore che si ottiene moltiplicando scalarmente
     il quadrato di omega[1] per l'autovettore v2 */
  for(i=0;i<2*n;i++)
    vtmp[i] += omega[1] * omega[1] * v2[i];
  printf("  Verifichiamo che v2 e' effettivamente un autovettore della matrice\n");
  printf("  JQJQ per l'autovalore -omega[1]*omega[1] a meno del seguente errore:\n");
  printf("    %le\n", norma_euclidea(vtmp, 2*n) );
  /* Definisco il vettore vtmp in modo tale che vtmp = J v2 */
  prod_matr_vet(J, v2, 2*n, vtmp);
  printf("  Verifichiamo che il prodotto scalare tra l'autovettore v1 e \n");
  printf("  il vettore J * v2 e' uguale a -1 a meno del seguente errore:\n");
  printf("    %le\n", fabs( prod_scal(v1, vtmp, 2*n) + 1 ) );
  /* Poniamo la seconda e quarta colonna di V uguali, rispettivamente
     a v1 e v2 */
  for(i=0;i<2*n;i++) {
    V[i][1] = v1[i];
    V[i][3] = v2[i];
  }
  /* Inizio del test che verifica che V e' simplettica */
  /* Eseguiamo il prodotto matriciale JV = J * V, facendo ricorso alla
     apposita function prod_matr */
  prod_matr(J, V, 2*n, JV);
  /* Poniamo S uguale alla matrice trasposta di V, facendo ricorso alla
     apposita function trasp_matr */
  trasp_matr(V, 2*n, S);
  /* Eseguiamo il prodotto matriciale VtJV = S * JV, facendo ricorso alla
     apposita function prod_matr */
  prod_matr(S, JV, 2*n, VtJV);
  /* Poniamo la matrice VtJV_men_J uguale alla matrice VtJV meno la
     matrice J (tale operazione matriciale deve essere svolta,
     SEPARATAMENTE, per ciascun elemento) */
  for(i=0; i<2*n; i++)
    for(j=0; j<2*n; j++)
      VtJV_men_J[i][j] = VtJV[i][j] - J[i][j];
  /* Inizialmente, azzeriamo il valore di max_VtJV_men_J */
  max_VtJV_men_J = 0;
  /* Poniamo il valore di max_VtJV_men_J uguale al massimo valore
     assoluto di TUTTI gli elementi di VtJV_men_J */
  for(i=0; i<2*n; i++)
    for(j=0; j<2*n; j++)
      if( max_VtJV_men_J < fabs(VtJV_men_J[i][j]) )
	max_VtJV_men_J = fabs(VtJV_men_J[i][j]);
  /* Stampiamo il massimo valore assoluto di tutti gli elementi di
     VtJV_men_J */
  printf("  La matrice V (ottenuta accostando gli autovettori precedentemente\n");
  printf("  calcolati della matrice JQJQ), e' simplettica a meno del seguente\n");
  printf("  errore massimo:\n    %le\n\n\n", max_VtJV_men_J);
  /* Fine del test che verifica che V e' simplettica */




  printf("\n\n\n     ###   Obiettivo 6   ###     \n");
  double QV[NDIM][NDIM], VtQV[NDIM][NDIM];
  double beta[NDIM/2], L[NDIM][NDIM], W[NDIM][NDIM];
  double QW[NDIM][NDIM], WtQW[NDIM][NDIM];

  /* Eseguiamo il prodotto matriciale QV = Q * V, facendo ricorso alla
     apposita function prod_matr */
  prod_matr(Q, V, 2*n, QV);
  /* Poniamo S uguale alla matrice trasposta di V, facendo ricorso alla
     apposita function trasp_matr */
  trasp_matr(V, 2*n, S);
  /* Eseguiamo il prodotto matriciale VtQV = S * QV, facendo ricorso alla
     apposita function prod_matr */
  prod_matr(S, QV, 2*n, VtQV);
  /* Calcolo delle componenti del vettore beta */
  for(i=0;i<n;i++)
    beta[i] = - VtQV[i][i+n] / VtQV[i][i];
  /* Definiamo la matrice L che consente di mettere VtQV in forma
     diagonale */
  for(i=0; i<2*n; i++)
    for(j=0; j<2*n; j++)
      if(i==j)
	L[i][j] = 1.;
      else
	L[i][j] = 0.;
  for(i=0; i<n; i++)
    L[i][i+n] = beta[i];
  /* Eseguiamo il prodotto matriciale W = V * L, facendo ricorso alla
     apposita function prod_matr */
  prod_matr(V, L, 2*n, W);
  /* Eseguiamo il prodotto matriciale QW = QW, facendo ricorso alla
     apposita function prod_matr */
  prod_matr(Q, W, 2*n, QW);
  /* Poniamo S uguale alla matrice trasposta di W, facendo ricorso alla
     apposita function trasp_matr */
  trasp_matr(W, 2*n, S);
  /* Eseguiamo il prodotto matriciale WtQW = S * QW, facendo ricorso alla
     apposita function prod_matr */
  prod_matr(S, QW, 2*n, WtQW);

  /* Stampa (ordinata) della matrice WtQW */
  printf("  Verifichiamo che la matrice W e' effettivamente in grado di\n");
  printf("  porre la matrice Q in forma diagonale. Infatti, il prodotto\n");
  printf("  matriciale tra la trasposta di W e QW e' uguale a:\n");
  for(i=0; i<2*n; i++) {
    printf("(");
    for(j=0; j<2*n; j++)
      printf(" %20.16lf", WtQW[i][j]);
    printf(" )\n");
  }




  printf("\n\n\n     ###   Obiettivo 7   ###     \n");
  double D[NDIM][NDIM], U[NDIM][NDIM], QU[NDIM][NDIM], UtQU[NDIM][NDIM];
  double JU[NDIM][NDIM], UtJU[NDIM][NDIM], UtJU_men_J[NDIM][NDIM];
  double max_UtJU_men_J;

  /* Definiamo la matrice D che consente di mettere WtQW in forma
     canonica */
  for(i=0; i<2*n; i++)
    for(j=0; j<2*n; j++)
      D[i][j] = 0.;
  for(i=0; i<n; i++) {
    D[i][i] = pow( WtQW[i+n][i+n] / WtQW[i][i] , 0.25 );
    D[i+n][i+n] = pow( WtQW[i][i] / WtQW[i+n][i+n] , 0.25 );
  }
  /* Eseguiamo il prodotto matriciale U = W * D, facendo ricorso alla
     apposita function prod_matr */
  prod_matr(W, D, 2*n, U);
  /* Eseguiamo il prodotto matriciale QU = QU, facendo ricorso alla
     apposita function prod_matr */
  prod_matr(Q, U, 2*n, QU);
  /* Poniamo S uguale alla matrice trasposta di U, facendo ricorso alla
     apposita function trasp_matr */
  trasp_matr(U, 2*n, S);
  /* Eseguiamo il prodotto matriciale UtQU = S * QU, facendo ricorso alla
     apposita function prod_matr */
  prod_matr(S, QU, 2*n, UtQU);

  /* Stampa (ordinata) della matrice UtQU */
  printf("  Verifichiamo che la matrice U e' effettivamente in grado di\n");
  printf("  porre la matrice Q in forma canonica (cioe', con gli autovalori\n");
  printf("  ripetuti due volte sulla diagonale e tutti gli altri elementi\n");
  printf("  quasi uguali a zero, a meno degli errori numerici). Infatti, il\n");
  printf("  prodotto matriciale tra la trasposta di U e QU e' uguale a:\n");
  for(i=0; i<2*n; i++) {
    printf("(");
    for(j=0; j<2*n; j++)
      printf(" %20.16lf", UtQU[i][j]);
    printf(" )\n");
  }
  /* Inizio del test che verifica che U e' simplettica */
  /* Eseguiamo il prodotto matriciale JU = J * U, facendo ricorso alla
     apposita function prod_matr */
  prod_matr(J, U, 2*n, JU);
  /* Poniamo S uguale alla matrice trasposta di U, facendo ricorso alla
     apposita function trasp_matr */
  trasp_matr(U, 2*n, S);
  /* Eseguiamo il prodotto matriciale UtJU = S * JU, facendo ricorso alla
     apposita function prod_matr */
  prod_matr(S, JU, 2*n, UtJU);
  /* Poniamo la matrice UtJU_men_J uguale alla matrice UtJU meno la
     matrice J (tale operazione matriciale deve essere svolta,
     SEPARATAMENTE, per ciascun elemento) */
  for(i=0; i<2*n; i++)
    for(j=0; j<2*n; j++)
      UtJU_men_J[i][j] = UtJU[i][j] - J[i][j];
  /* Inizialmente, azzeriamo il valore di max_UtJU_men_J */
  max_UtJU_men_J = 0;
  /* Poniamo il valore di max_UtJU_men_J uguale al massimo valore
     assoluto di TUTTI gli elementi di UtJU_men_J */
  for(i=0; i<2*n; i++)
    for(j=0; j<2*n; j++)
      if( max_UtJU_men_J < fabs(UtJU_men_J[i][j]) )
	max_UtJU_men_J = fabs(UtJU_men_J[i][j]);
  /* Stampiamo il massimo valore assoluto di tutti gli elementi di
     UtJU_men_J */
  printf("  La matrice U e' simplettica a meno del seguente ");
  printf("errore massimo:\n    %le\n\n\n", max_UtJU_men_J);
  /* Fine del test che verifica che U e' simplettica */



}


/* ### */
double determinante(double mat[][NDIM], int n)
{
  /* QUESTA FUNCTION E' ESATTAMENTE IDENTICA A QUELLA CHE SI TROVA (ad
     esempio) IN prod_det.c */
  int i;
  double ris=0., mattmp[NDIM][NDIM];
  /* Questa funzione effettua il calcolo del determinante di una
     matrice nxn-dimensionale utilizzando il metodo (ricorsivo)
     di Laplace */
  /* Se n e' maggiore di 1 ... */
  if(n > 1) {
    /* ... Facciamo un ciclo per i che va da 0 a n-1 */
    for(i=0; i<n; i++) {
      /* Mettiamo nella matrice mattmp una parte della matrice mat,
	 che otteniamo da mat stessa togliendo la 0-esima riga e la
	 i-esima colonna. Facciamo tutto cio', chiamando la function
	 togli_riga0_colonnai */
      togli_riga0_colonnai(mat, n, i, mattmp);
      /* Aggiorniamo il risultato, sommandovi il prodotto tra:
	 (1) l'elemento della matrice mat che sta nella 0-esima riga
	     e nella i-esima colonna,
	 (2) il minore corrispondente, ovvero il determinante della
	     matrice (n-1)x(n-1)-dimensionale che e' posta nella matrice
	     mattmp,
	 (3) il segno che si ottiene elevando -1 alla i    */	 
      ris += ( 1 - 2 * (i % 2) ) * mat[0][i] * determinante(mattmp, n-1);
    }
  }
  /* Altrimenti ... */
  else {
    /* ... Poniamo il risultato uguale al termine che sta nel vertice
       sinistro in alto della matrice mat */
    ris = mat[0][0];
  }
  /* Restituiamo il valore del risultato alla funzione chiamante */
  return(ris);
}


/* ### */
void gramschmidt_simpl_duevet(double v1[], double v2[], int n)
{
  int i, j;
  double vtmp[NDIM], J[NDIM][NDIM], norma, coef_prod_simpl;

  /* Definiamo J uguale alla matrice fondamentale della geometria
     simplettica */
  for(i=0; i<n; i++)
    for(j=0; j<n; j++)
      J[i][j] = 0.;
  for(i=0; i<n/2; i++)
    J[i][i+n/2] = -1;
  for(i=0; i<n/2; i++)
    J[i+n/2][i] = 1;

  /* Pongo norma uguale a || v1 ||, cioe' la norma euclidea del
     vettore v1 */
  norma = norma_euclidea(v1, n);
  /* Normalizzo il vettore v1 */
  for(i=0;i<n;i++)
    v1[i] /= norma;
  /* Definisco il vettore vtmp in modo tale che vtmp = J v2 */
  prod_matr_vet(J, v2, n, vtmp);
  /* Pongo la variabile coef_prod_simpl uguale al prodotto scalare tra
     i vettori v1 e vtmp.  Questo calcolo viene effettuato facendo
     ricorso alla apposita function prod_scal */
  coef_prod_simpl = prod_scal(v1, vtmp, n);
  /* Divido ogni i-esimo elemento di v2 per il valore di
     coef_prod_simpl cambiato di segno */
  for(i=0; i<n; i++)
    v2[i] /= -coef_prod_simpl;

}


/* ### */
void inv_mat_cramer(double A[][NDIM], int n, double A_inv[][NDIM])
{
  int i, j, l, m;
  double detA, b[NDIM], mattmp[NDIM][NDIM];
  /* Questa funzione calcola l'inversa della matrice A, utilizzando il
     metodo di Cramer  */
  /* Calcolo del determinante della matrice A, facendo ricorso alla
     function determinante */
  detA = determinante(A, n);
  /* Eseguo un ciclo per i che va da 0 a n-1 */
  for(i=0; i<n; i++) {
    /* Azzero la i-esima componente del "vettore noto" b */
    b[i] = 0.;
  }
  /* Eseguo un ciclo per m che va da 0 a n-1 */
  for(m=0; m<n; m++) {
    /* Pongo la m-esima componente del "vettore noto" b uguale a 1 */
    b[m] = 1.;
    /* Eseguo un ciclo per i che va da 0 a n-1 */
    for(i=0; i<n; i++) {
      /* Eseguo un ciclo per j che va da 0 a n-1 */
      for(j=0; j<n; j++) {
	/* Eseguo un ciclo per l che va da 0 a n-1 */
	for(l=0; l<n; l++) {
	  /* Se i e' uguale a l, poniamo l'elemento che occupa la j-esima
	     riga e la l-esima colonna della matrice mattmp uguale alla
	     j-esima componente del "vettore noto" b ... */
	  if(i == l)
	    mattmp[j][l] = b[j];
	  /* ...Altrimenti, poniamo l'elemento che occupa la j-esima
	     riga e la l-esima colonna della matrice mattmp uguale al
	     corrispondente elemento della matrice A */
	  else
	    mattmp[j][l] = A[j][l];
	}
      }
      /* Poniamo l'elemento che occupa la i-esima riga e la m-esima
	 colonna della matrice A_inv uguale al determinante della
	 matrice mattmp diviso per il determinante di A */
      A_inv[i][m] = determinante(mattmp, n) / detA;
    }
    /* Azzero la m-esima componente del "vettore noto" b */
    b[m] = 0.;
  }
}


/* ### */
double met_potenze_simpl(double A[][NDIM], int n, double v1[], double v2[])
{
  int i,j;
  double w1[NDIM], w2[NDIM];
  double lambda, lambdaold, norma;

  /* Definisco l'approssimazione iniziale dell'autovalore lambda in
     modo tale che sia uguale a 0 */
  lambda = 0.;
  /* Definisco l'approssimazione iniziale dell'autovettore v1 in modo
     tale che v1 abbia norma uguale a 1 e tutte le sue componenti siano
     uguali tra loro */
  v1[0] = 1. / sqrt(n);
  for(i=1; i<n; i++)
    v1[i] = v1[0];
  /* Definisco l'approssimazione iniziale dell'autovettore v2 in modo
     tale che v2 abbia tutte le sue componenti uguali a quelle di v1,
     fatta eccezione per la prima componente che avra' segno
     opposto */
  for(i=0; i<n; i++)
    v2[i] = v1[i];
  v2[0] = -v1[0];

  /* Invoco la funzione che effettua l'"ortonormalizzazione
     simplettica" alla Gram-Schmidt per i due soli vettori v1 e v2 */
  gramschmidt_simpl_duevet(v1, v2, n);
  /* Eseguiamo il ciclo di operazioni seguenti mentre e' verificata
     la condizione che scriviamo dopo la prossima parentesi graffa
     chiusa (cioe' } ) */
  do {
    /* Ridefinisco la vecchia approssimazione dell'autovalore uguale 
       a quella attuale */
    lambdaold = lambda;
    /* Calcolo del vettore w1 uguale al prodotto di A per v1,
       facendo ricorso alla function prod_matr_vet */
    prod_matr_vet(A, v1, n, w1);
    /* La nuova approssimazione dell'autovalore lambda e' data dal
       prodotto scalare v1 per w1 */
    lambda = prod_scal(v1, w1, n);
    /* Ridefiniamo v1 in modo che sia uguale a w1 */
    for(i=0;i<n;i++)
      v1[i] = w1[i];
    /* Calcolo del vettore w2 uguale al prodotto di A per v2,
       facendo ricorso alla function prod_matr_vet */
    prod_matr_vet(A, v2, n, w2);
    /* Ridefiniamo v2 in modo che sia uguale a w2 */
    for(i=0;i<n;i++)
      v2[i] = w2[i];
    /* Invoco la funzione che effettua l'"ortonormalizzazione
       simplettica" alla Gram-Schmidt per i due soli vettori v1 e v2 */
    gramschmidt_simpl_duevet(v1, v2, n);
    /* Dopo la parentesi graffa chiusa seguente eseguiamo il test che
       controlla l'iterazione del ciclo: se il rapporto tra il valore
       assoluto | lambdaold - lambda | e il valore assoluto di lambda
       e' superiore a 2n * 10 elevato alla -16 (cioe' circa n volte
       l'errore di macchina), continuiamo ad iterare l'esecuzione del
       ciclo. */
  } while( fabs(lambdaold -lambda) / fabs( lambda ) > n * 2.e-16 );

  return(lambda);
}


/* ### */
double norma_euclidea(double v[], int n)
{
  int i;
  double norma = 0;
  for(i=0; i<n; i++)
    norma += v[i] * v[i];
  norma = sqrt(norma);
  return(norma);
}


/* ### */
void prod_matr(A, B, n, C)
     double A[][NDIM], B[][NDIM], C[][NDIM];
     int n;
{
  int i, j;
  double Bt[NDIM][NDIM];

  /* Pongo Bt uguale alla matrice trasposta di B, facendo ricorso alla
     apposita function trasp_matr */
  trasp_matr(B, n, Bt);

  for(i=0; i<n; i++)
    for(j=0; j<n; j++)
      /* Pongo l'elemento ij-esimo di C uguale al prodotto scalare tra
	 la i-esima riga e la j-esima riga di Bt (che e' uguale alla
	 j-esima colonna di B). Questo calcolo viene effettuato
	 facendo ricorso alla apposita function prod_scal */
      C[i][j] = prod_scal(&A[i][0] , &Bt[j][0], n);

}


/* ### */
void prod_matr_vet(double mat[][NDIM],double vet[], int n, double vetris[])
{
  int i, j;
  /* Questa funzione effettua il calcolo riga per colonna tra la
     matrice mat (nxn-dimensionale) e il vettore vet (n-dimensionale),
     il risultato e' posto nella matrice vetris  */
  for(i=0;i<n;i++) {
    vetris[i] = 0.;
    for(j=0;j<n;j++) {
      vetris[i] += mat[i][j] * vet[j];
    }
  }
}


/* ### */
double prod_scal(v1, v2, n)
     double v1[], v2[];
     int n;
{
  int i;
  double ris = 0.;
  for(i=0; i<n; i++)
    ris += v1[i] * v2[i];
  return(ris);
}


/* ### */
int risolvi_eq(a, b, c, x1, x2)
     double a, b, c, *x1, *x2;
{
  int flag;
  double Delta;
  /* Poniamo Delta = b^2 - 4*a*c */
  Delta = b*b - 4.*a*c;
  /* Distinguiamo i diversi casi a seconda del segno di Delta */
  if(Delta > 0.) {
    /* Siamo nel caso delle due radici reali e distinte:
       calcoliamo x1 ... */
    *x1 = -0.5 * ( b + sqrt(Delta) ) / a;
    /* ... e x2 */
    *x2 = -0.5 * ( b - sqrt(Delta) ) / a;
    /* Ricordiamoci di inizializzare flag */
    flag = 1;
  }
  else if(Delta == 0.) {
    /* Siamo nel caso delle due radici reali coincidenti:
       calcoliamo x1 e x2 */
    *x1 = *x2 = -0.5 * b / a;
    /* Ricordiamoci di inizializzare flag */
    flag = 0;
  }
  else {
    /* Siamo nel caso delle due radici complesse:
       poniamo x1 uguale alla parte reale delle soluzioni ... */
    *x1 = -0.5 * b / a;
    /* ... e poniamo x2 uguale alla parte immaginaria delle soluzioni */
    *x2 = 0.5 * sqrt(-Delta) / a;
    /* Ricordiamoci di inizializzare flag */
    flag = -1;
  }
  /* Usciamo dalla funzione restituendo il valore di flag */
  return(flag);
}


/* ### */
void togli_riga0_colonnai(mat, n, i, matridotta)
     double mat[][NDIM], matridotta[][NDIM];
     int n, i;
{
  /* QUESTA FUNCTION E' ESATTAMENTE IDENTICA A QUELLA CHE SI TROVA (ad
     esempio) IN prod_det.c */
  int j, l;
  /* Questa funzione isola una parte della matrice mat, che otteniamo
     da mat stessa togliendo la 0-esima riga e la i-esima colonna, il
     risultato e' posto nella matrice matridotta; mat e' da
     considerarsi nxn-dimensionale e matridotta
     (n-1)x(n-1)-dimensionale */
  for(j=1;j<n;j++) {
    for(l=0;l<n;l++) {
      if(l!=i)
	matridotta[j-1][l-(l>i)]=mat[j][l];
    }
  }
}


/* ### */
void togli_rigal_colonnal(mat, n, l, matridotta)
     double mat[][NDIM], matridotta[][NDIM];
     int n, l;
{
  /* QUESTA FUNCTION E' UNA LIEVISSIMA VARIAZIONE DELLA FUNCTION
     PRECEDENTE, cioe' togli_riga0_colonnai */
  int i, j;
  /* Questa funzione isola una parte della matrice mat, che otteniamo
     da mat stessa togliendo la l-esima riga e la l-esima colonna, il
     risultato e' posto nella matrice matridotta; mat e' da
     considerarsi nxn-dimensionale e matridotta
     (n-1)x(n-1)-dimensionale */
  for(i=0;i<n;i++) {
    for(j=0;j<n;j++) {
      if(i!=l && j!=l)
	matridotta[i-(i>l)][j-(j>l)]=mat[i][j];
    }
  }
}


/* ### */
void trasp_matr(A, n, At)
     double A[][NDIM], At[][NDIM];
     int n;
{
  int i, j;

  /* Pongo At uguale alla matrice trasposta di A (tale operazione
     matriciale deve essere svolta, SEPARATAMENTE, per ciascun
     elemento) */
  for(i=0; i<n; i++)
    for(j=0; j<n; j++)
      At[i][j] = A[j][i];

}