// funzione f(x) da interpolare
// nell'intervallo a,b
double a=-5,b=5;
double f(double x) {
  //return x*sin(x);
  return 1./(1.+x*x);
};

// Numero di nodi da usare nell'interpolazione
   int N=16;
// passo per interpolazione (punti equidistanti)
   double h=(b-a)/(N-1.);
// vettori che conterranno i nodi dell'intepolazione
   double x[100],y[100];
// funzione interpolante.
// Solo definita l'implementazione è in fondo
// al listato
double g(double X);

//main
void Interpolazione2() {
    gStyle->SetOptStat(0);
// calcolo dei nodi
   printf("%f\n",h);
   for (int i=0; i<N ; i++) {
       x[i]=a+h*i;
       y[i]=f(x[i]);
       printf("%f %f\n",x[i],y[i]);

   }
// grafico della funzione "vera" per confronto
   //TH2F* hb = new TH2F("hbox","",100,-6.,6.,100,-6.,2.);
    TH2F* hb = new TH2F("hbox","",100,-6.,6.,100,-0.5,2.5);
    hb->Draw();
   TF1 * f=new TF1("f","f(x)",a,b);
   f->SetLineWidth(5);
   f->Draw("same");
    //f->Draw();
// Grafico della funzione interpolante
   TF1 * g=new TF1("g","g(x)",a,b);
   g->SetLineColor(kBlue);
   g->SetLineWidth(4.0);
   g->SetLineStyle(1);
   g->Draw("same");
// grafico dei punti usati per l'interpolazione
   TGraph * p=new TGraph(N,x,y);
   p->SetMarkerStyle(21);
   p->SetMarkerSize(0.8);
   p->SetMarkerColor(1);
   p->Draw("Psame");
   for (int i=0; i<N ; i++) {
       x[i]=a+h*i;
       double yy=g->Eval(x[i]);
       printf("Funzione interpolante a x=%f vale %f errore al nodo = %f\n",x[i],yy, yy-y[i]);
   }

}
double g(double X){
// interpolazione con polinomio di grado N nel punto x
// converto x in t
   double t=(X-a)/h;
// vettore contenente il valore dei polinomi di Lagrange
   double L[100];
   for (int i=0;i<N;i++) {
// numeratore polinomio i-esimo
       double num=1.;
// denonominatore polinomio i-esimo
       double den=1.;
       for (int j=0; j<N; j++) {
          // escludo caso j=i
          if (j==i) continue;
          num=num*(t-j);
          den=den*(i-j);
       }
       L[i]=num/den;
   }
// calcolo del valore della funzione interpolante nel punto x
   double g=0;
   for (int i=0; i<N;i++) {
       g=g+y[i]*L[i];
   }
   return g;
}