/* pol_legendre.c : file sorgente di un programma che consente di studiare
                    alcune caratteristiche dei polinomi di Legendre */

#include<stdio.h>
#include<math.h>
/* E' necessario includere il file stdlib.h, al fine di poter
   utilizzare la function exit */
#include<stdlib.h>

/* Definizione di alcuni parametri utili */
#define NMAX 200

/* Definizione dei prototipi delle funzioni */
void calc_pol_leg(int n, double P_n[]);
void copia_polinomio(double coef_pol_1[], int grado_pol, double coef_pol_2[]);
void deriv_polinomio(double coef_pol[], int grado_pol, double coef_der_pol[]);
void input_polinomio(double coef_pol[], int grado_pol);
void prod_pol_per_scal(double coef_pol[], int grado_pol, double coef_scal);
void prod_pol_per_x(double coef_pol[], int grado_pol);
void somma_polinomi(double coef_pol_add_1[], int grado_pol_1,
		    double coef_pol_add_2[], int grado_pol_2,
		    double coef_pol_ris[]);
void stampa_polinomio(double coef_pol[], int grado_pol);
double valuta_polinomio(double coef_pol[], int grado_pol, double x);

int main()
{


  printf("\n\n     ###   Obiettivo 1   ###     \n");
  int n;
  double P_n[NMAX+1], P_nmen1[NMAX], P_nmen2[NMAX-1], x, err_test;

  do {
    printf("  Valore di n per cui si vuole verificare l'identita' ");
    printf("ricorsiva? (2<=n<=%d)\n  ", NMAX);
    scanf("%d", &n);
  } while(n<2 || n > NMAX);

  printf("\n     *** Input del polinomio P_%d(x) ***\n", n-2);
  input_polinomio(P_nmen2, n-2);
  printf("\n     *** Input del polinomio P_%d(x) ***\n", n-1);
  input_polinomio(P_nmen1, n-1);
  printf("\n     *** Input del polinomio P_%d(x) ***\n", n);
  input_polinomio(P_n, n);

  printf("\n  Inserisci un valore di x :  ");
  scanf("%lf", &x);
  printf("\n  Risultato di %d x P_%d(x) - %d P_%d(x) - %d P_%d(x) ",
	 2*n-1, n-1, n-1, n-2, n, n);
  printf("quando x = %lf :\n", x);
  err_test = (2*n-1) * x * valuta_polinomio(P_nmen1,n-1,x) -
             (n-1) * valuta_polinomio(P_nmen2,n-2,x) -
             n * valuta_polinomio(P_n,n,x);
  printf("    %le\n", err_test);


  printf("\n\n     ###   Obiettivo 2   ###     \n");
  int esp_10;
  double t;

  if(n != 2 || fabs(x) > 1) {
    printf("\n  La verifica richiesta dall'obiettivo 2 non puo' essere\n");
    printf("  effettuata, perche' devono essere verificate ENTRAMBE\n");
    printf("  le seguenti condizioni:\n");
    printf("      n = 2     e     |x| <= 1\n");
  }
  else {
    printf("\n  Verifica empirica della seguente relazione : \n");
    printf("  P_0(x) + P_1(x) t + P_2(x) t^2 - 1 / sqrt( 1 - 2 x t + t^2 ) ");
    printf("= O( t^3 )\n");
    printf("  -------------------------------------------------------------");
    printf("----------\n");
    printf("        t                errore\n");
    printf("  __________________________________\n");
    for(esp_10 = 0; esp_10 <= 5; esp_10++) {
      t = 0.9 / pow(10, esp_10);
      err_test = valuta_polinomio(P_nmen2,0,x) +
	         valuta_polinomio(P_nmen1,1,x) * t +
	         valuta_polinomio(P_n,2,x) * t * t;
      err_test -= 1 / sqrt( 1 - 2*x*t + t*t );
      err_test = fabs(err_test);
      printf("    %le      %le\n", t, err_test);
    }
  }


  printf("\n\n     ###   Obiettivo 3   ###     \n");
  double prima_deriv[NMAX], seconda_deriv[NMAX+1], prod_pol_der[NMAX+2];
  double pol_tmp[NMAX+1], S[NMAX+1];

  deriv_polinomio(P_n, n, prima_deriv);
  copia_polinomio(prima_deriv, n-1, prod_pol_der);
  prod_pol_per_x(prod_pol_der, n-1);
  prod_pol_per_x(prod_pol_der, n);
  prod_pol_per_scal(prod_pol_der, n+1, -1);
  somma_polinomi(prod_pol_der, n+1, prima_deriv, n-1, pol_tmp);
  deriv_polinomio(pol_tmp, n+1, seconda_deriv);
  copia_polinomio(P_n, n, pol_tmp);
  prod_pol_per_scal(pol_tmp, n, n*(n+1));
  somma_polinomi(seconda_deriv, n, pol_tmp, n, S);

  printf("\n  Risultato di d/dx [ (1 - x^2) d/dx P_%d(x) ] - %d P_%d(x) ",
	 n, n*(n+1), n);
  printf("(per x generico) : \n");
  stampa_polinomio(S, n);


  printf("\n\n     ###   Obiettivo 4   ###     \n");
  do {
    printf("  Valore di n per cui si vuole calcolare l'espansione del\n");
    printf("  corrispondente polinomio di Legendre P_n(x)? (2<=n<=%d)\n  ",
	   NMAX);
    scanf("%d", &n);
  } while(n<2 || n > NMAX);

  calc_pol_leg(n, P_n);
  printf("\n  Espansione del polinomio P_%d(x) :\n", n);
  stampa_polinomio(P_n, n);


  printf("\n\n     ###   Obiettivo 5   ###     \n");
  double t_alla_n;

  do {
    printf("  Inserire un valore di x (-1 <= x <= 1) :  ");
    scanf("%lf", &x);
  } while(x<-1 || x>1);
  do {
    printf("  Inserire un valore di t (0 <= t <= 1/2) :  ");
    scanf("%lf", &t);
  } while(t<0 || t>0.5);

  err_test = 0;
  t_alla_n = 1;
  for(n=0;n<=NMAX;n++) {
    calc_pol_leg(n, P_n);
    err_test += valuta_polinomio(P_n,n,x) * t_alla_n;
    t_alla_n *= t;
  }
  err_test -= 1 / sqrt( 1 - 2*x*t + t*t );

  printf("\n  Verifica dell'espansione in serie di ");
  printf("1 / sqrt( 1 - 2 x t + t^2 ) : \n");
  printf("    %le\n\n\n", err_test);


}


/* ### */
void calc_pol_leg(n, P_n)
  int n;
  double P_n[];
{
  int i;
  double P_i[NMAX+1], P_imen1[NMAX], P_imen2[NMAX-1];
  double pol_tmp_1[NMAX+1], pol_tmp_2[NMAX+1];

  P_imen2[0] = 1;
  P_imen1[0] = 0;  P_imen1[1] = 1;

  for(i=2;i<=n;i++) {
    copia_polinomio(P_imen1, i-1, pol_tmp_1);
    prod_pol_per_scal(pol_tmp_1, i-1, 2*i-1);
    prod_pol_per_x(pol_tmp_1, i-1);
    copia_polinomio(P_imen2, i-2, pol_tmp_2);
    prod_pol_per_scal(pol_tmp_2, i-2, -(i-1));
    somma_polinomi(pol_tmp_1, i, pol_tmp_2, i-2, P_i);
    prod_pol_per_scal(P_i, i, 1/((double) i));
    copia_polinomio(P_imen1, i-1, P_imen2);
    copia_polinomio(P_i, i, P_imen1);
  }

  if(n==0)
    copia_polinomio(P_imen2, 0, P_n);
  else if(n==1)
    copia_polinomio(P_imen1, 1, P_n);
  else
    copia_polinomio(P_i, n, P_n);
  return;
}

/* ### */
void copia_polinomio(coef_pol_1, grado_pol, coef_pol_2)
  double coef_pol_1[];
  int grado_pol;
  double coef_pol_2[];
{
  int i;
  for(i=0; i<=grado_pol; i++)
    coef_pol_2[i] = coef_pol_1[i];
}

/* ### */
void deriv_polinomio(coef_pol, grado_pol, coef_der_pol)
     double coef_pol[];
     int grado_pol;
     double coef_der_pol[];
{
  int i;
  for(i=1; i<=grado_pol; i++)
    coef_der_pol[i-1] = i * coef_pol[i];
}

/* ### */
void input_polinomio(coef_pol, grado_pol)
     double coef_pol[];
     int grado_pol;
{
  int i, denom;
  printf("  Denominatore?  ");
  scanf("%d", &denom);
  for(i=grado_pol;i>=0;i--) {
    if( i%2 == grado_pol%2) {
      printf("  Coef. di x elevato alla %2d :  ", i);
      scanf("%lf", coef_pol+i);
      coef_pol[i] /= denom;
    }
    else
      coef_pol[i] = 0;
  }
}

/* ### */
void prod_pol_per_scal(coef_pol, grado_pol, coef_scal)
     double coef_pol[];
     int grado_pol;
     double coef_scal;
{
  int i;
  for(i=0;i<=grado_pol;i++)
    coef_pol[i] *= coef_scal;
  return;
}

/* ### */
void prod_pol_per_x(coef_pol, grado_pol)
     double coef_pol[];
     int grado_pol;
{
  int i;
  for(i=grado_pol;i>=0;i--)
    coef_pol[i+1] = coef_pol[i];
  coef_pol[0] = 0;
  return;
}

/* ### */
void somma_polinomi(coef_pol_add_1, grado_pol_1, coef_pol_add_2,
		    grado_pol_2, coef_pol_ris)
     double coef_pol_add_1[];
     int grado_pol_1;
     double coef_pol_add_2[];
     int grado_pol_2;
     double coef_pol_ris[];
{
  int i;
  if(grado_pol_1 < grado_pol_2) {
    printf("  La funzione somma_polinomi e' stata scritta in modo tale da\n");
    printf("  poter funzionare correttamente solo quando il grado del primo\n");
    printf("  polinomio addendo e' superiore o uguale al grado del secondo!\n");
    exit(1);
  }
  for(i=0;i<=grado_pol_2;i++)
    coef_pol_ris[i] = coef_pol_add_1[i] + coef_pol_add_2[i];
  for(i=grado_pol_2+1;i<=grado_pol_1;i++)
    coef_pol_ris[i] = coef_pol_add_1[i];
  return;
}

/* ### */
void stampa_polinomio(coef_pol, grado_pol)
     double coef_pol[];
     int grado_pol;
{
  int i;
  for(i=0;i<=grado_pol;i++)
    printf("  Coef. di x elevato alla %2d :  %le\n", i, coef_pol[i]);
  return;
}


/* ### */
double valuta_polinomio(coef_pol, grado_pol, x)
     double coef_pol[];
     int grado_pol;
     double x;
{
  int i;
  double potx=1, ris=0.;
  for(i=0;i<=grado_pol;i++) {
    ris += coef_pol[i] * potx;
    potx *= x;
  }
  return ris;
}