// funzione f(x) da interpolare
// nell'intervallo a,b
double a=-5,b=5;
double f(double x) {
  //   return 1./(1.+x*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, bool verbose=false);
double gComposto(double X, int Nnodi, bool verbose=false);
void Interpolazione5() {
// 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]);

   }

   gStyle->SetOptStat(0);
// Grafico della funzione interpolante
   std::string fname="gComposto(x,16,false)";
   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 * ftrue=new TF1("f","f(x)",a,b);
   ftrue->Draw("same");
   //return;

// analizziamo una posizione estrema
   //double X=b-h/2.;
// analizziamo una posizione vicina a 0 
   //double X=0.235;
   double X=0.85;
   double gC = gComposto(X,N,true);
   std::cout<<"gComposto ("<<X<<") = "<< gC <<" errore = "<<gC-f(X)<<std::endl;
   //return;

   TF1 * g=new TF1("g",fname.c_str(),a,b);
   /*
   std::cout<<" Loop - 100 punti - gComposto ..... "<<std::endl;
   for (int i=0; i<100; ++i)
     {
       X= a+i*(b-a)/100+0.005;
       if (X > b) continue;
       gC = gComposto(X,N,true);
       std::cout<<i<<"  gComposto("<<X<<") = "<<g->Eval(X) <<std::endl;
     }
   //return;
   */
   g->SetLineColor(4);
   g->SetLineWidth(4.0);
   g->Draw("same");
   
   

}
double gComposto(double X, int Nnodi, bool verbose=false)
{
   double ginterp=-999.;
   double ginterpprev = -999.;
   double diffip=9999999.;
   double diffi =9999999.;
   int N = Nnodi;
   if (verbose) std::cout<<"------------ calcolo gComposto(x = "<<X<<") con Nnodi = "<<N<<" in totale "<<std::endl;
   for (int j=1; j<N; j++) {
     ginterp = g(X,j,verbose);
     diffi=ginterp-ginterpprev;
     if (verbose) printf("iterazione con #Nodi = %d g(%f)= %f diff w.r.t iter. precedente = %f ; g, diff precedenti: %f, %f\n ",
			 j,X,ginterp,diffi, ginterpprev, diffip);
     if (j>2)
       {
	 if (fabs(diffi) > fabs(diffip) 
	     //     || double(j)>0.75*double(N)  
	     //     || (fabs(diffi)/fabs(ginterp)<0.00001 && j>2)
	     ) {
	   if (verbose) printf("@ #Nodi = %d g(%f)= %f diff w.r.t previous iteration = %f .... stop: migliore interp = %f con pol di grado %d\n ",j,X,ginterp,ginterp-ginterpprev, ginterpprev, j-2);
	   return ginterpprev;
	   break;
	 }
       }
     ginterpprev = ginterp;
     //diffip = diffi;
     //if (fabs(diffi/ginterp)>0.001 && j%2==0) diffip = diffi;
     if (fabs(diffi/ginterp)>0.001) diffip = diffi;
   }
   if (verbose) printf("No more nodes\n");
   return ginterp;
}
double g(double X, int grado, bool verbose=false){
// 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) {
     //t piu' vicino a i2
       i3=i2;
       i1=i3-(grado+1)/2;
   }
   else {
       i3=i1;
       i1=i3-grado/2;
   }

// fuori range se i1<0 o i2>N-1    
// 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;
   }

   if (verbose)   std::cout<<" i1, i2 = "<<i1<<" "<<i2<<" t="<<t<< std::endl; 
// 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];
       if (verbose) std::cout<<" L[i] = "<<L[i]<<" y[i]="<<y[i]<<" g = "<<g<<std::endl;
   }
   return g;
}