// funzione f(x) da interpolare
// nell'intervallo a,b
double a=-5,b=5;
double f(double 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.);
// vettore 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, int grado);
void Interpolazione4() {
// calcolo dei nodi
   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 interpolante
   std::string fname="gComposto(x,16)";
   TH2F* hbox = new TH2F("hbox",fname.c_str(),100,-6.,6.,100,0.,1.05);
   hbox->Draw();
   // 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");
// grafico della funzione "vera" per confronto
   TF1 * f=new TF1("f","f(x)",a,b);
   f->Draw("same");

// analizziamo una posizione estrema
   double X=b-h/2.;
   //double X=0.1;
   double ginterp=-999.;
   double ginterpprev = -999.;
   double diffip=9999999.;
   double diffi =9999999.;
   std::cout<<"------------ x = "<<X<<std::endl;
   for (int i=0; i<N-1; i++) {
     ginterp = g(X,i);
     if (i>0)
       {
	 diffi=ginterp-ginterpprev;
	 if (abs(diffi) > abs(diffip) || i>0.75*N) {
	   //printf("grado=%d g(x)=%f diff w.r.t previous iteration = %f .... stop: migliore interp = %f con pol di grado %d\n ",i,ginterp,ginterp-ginterpprev, ginterpprev, i-1);
	   //break;
	 }
       }
     printf("grado=%d g(x)=%f diff w.r.t previous iteration = %f\n",i,ginterp,ginterp-ginterpprev);
     ginterpprev = ginterp;
     if (abs(diffi)>0.00001) diffip = diffi;
   }

}
double g(double X, int grado){
// interpolazione con polinomio di grado N nel punto x
// converto x in t
   double t=(X-a)/h;
// identifico il punto più vicino
// t sarà compreso tra i1 e i2
   int i1,i2,i3;
   i1=floor(t);
   i2=i1+1;
// cerco chi dei due è il più vicino
// e identifico l'estremo inferiore del range
// di punti in base al grado
   if (t-i1>i2-t) {
       i3=i2;
       i1=i3-(grado+1)/2;
   }
   else {
       i3=i1;
       i1=i3-grado/2;
   }

// controllo se i1 è fuori range
   if (i1<0) i1=0;
// fisso i2 in base al grado
   i2=i1+grado;
// controllo se i2 è fuori range
   if (i2>N-1){
       i2=N-1;
       i1=i2-grado;
   }

// interpolazione di Lagrange locale
// vettore contenente il valore dei polinomi di Lagrance
   double L[100];
   for (int i=i1;i<=i2;i++) {
// numeratore polinomio i-esimo
       double num=1.;
// denonominatore polinomio i-esimo
       double den=1.;
       for (int j=i1; j<=i2; 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=i1; i<=i2;i++) {
       g=g+y[i]*L[i];
   }
   return g;
}