// 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 Interpolazione3() {
// 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);
   std::string fname="g(x,3)";
// Grafico della funzione interpolante
   TH2F* hbox = new TH2F("hbox",fname.c_str(),100,-6.,6.,100,0.,1.05);
   hbox->Draw();
   TF1 * g=new TF1("g",fname.c_str(),a,b);
   g->SetLineColor(4);
   g->SetLineWidth(4.0);
   g->Draw("same");

   std::string fname1="g(x,9)";
   TF1 * g5=new TF1("g",fname1.c_str(),a,b);
   g5->SetLineColor(4);
   g5->SetLineWidth(4.0);
   g5->SetLineStyle(2);
   g5->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");
// grafico della funzione "vera" per confronto
   TF1 * f=new TF1("f","f(x)",a,b);
   f->Draw("same");

   TLegend* tl = new TLegend(0.15,0.7,0.45,0.9);//xmin, ymin, xmax, ymax (% della pad)
  tl->SetBorderSize(0);
  tl->SetFillColor(0);
  tl->SetFillStyle(0);
  tl->SetTextFont(42);
  tl->SetTextSize(0.);
  tl->AddEntry(g,fname.c_str(),"L");
  tl->AddEntry(g5,fname1.c_str(),"L");
  tl->AddEntry(p,"Nodi","p");
  tl->AddEntry(f,"Funzione","L");
  tl->Draw();

   
}
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;
   }
   //std::cout<<"x= "<<X<<" i1, i2 "<<i1<<" "<<i2<<" NnodiLoc= "<<i2-i1+1<<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];
   }
   return g;
}