/*no se usa en el código pero una implementación a futuro podría ser implementar este tipo de interpolación
#include <cmath>
#include <iostream>
#include<iomanip>
//#include "algebralineal.h"
//#include "interpolacion.h"
//#include "gnuplot_i_calnum.h"
using namespace std;
using namespace CALNUM;

inline double fun(double x, int ft)
{
    switch(ft)
    {
    case 0:
        return 1./(1.+25*x*x);
        break;
    case 1:
        return abs(x);
        break;
    case 2:
        return sin(x);
    case 3:
        return 2*atan(x) +exp(5-x) + 5*cos(x);
    }
}

template <class T> int
splines(Vector<T> &x, Vector<T> &a, Vector<T> &b, Vector<T> &c, Vector<T> &d)

{
    cout<<"splines naturales"<<endl;

    int n=x.size(); //n numero de puntos
    int nint= n-1; //numero de intervalos
    int neq =nint -1; //numero de ecuaciones
    //Creacion vector diagonal superior
    Vector<T> h(nint);
    for( int i=0; i<nint; i++) h[i]=x[i+1]-x[i];
    //Creacion Vector diagonal
    Vector<T> u(nint), v(nint); //u(0), v(0) no estan definidos
    for( int i=1; i<nint; i++) u[i]=2*(h[i-1]+h[i]);
    for( int i=0; i<neq; i++)
        v[i+1]=3*( (a[i+2]-a[i+1])/h[i+1] - (a[i+1]-a[i])/h[i] );
    //Eliminacion de Gauss
    for( int i=1; i<neq; i++)
    {
        u[i+1]=u[i+1]-h[i]*h[i]/u[i];
        v[i+1]=v[i+1]-h[i]*v[i]/u[i];
    }
    //Sustitucion hacia atras
    c[nint-1]=v[nint-1]/u[nint-1];
    d[nint-1]= -c[nint-1]/(3*h[nint-1]);
    for( int k=nint-2; k>0; k--)
    {
        c[k]=(v[k]-h[k]*c[k+1])/u[k];
        d[k]=(c[k+1]-c[k])/(3*h[k]);
    }
    //Valores especiales en x0
    c[0]= T(0);
    d[0]=c[1]/(3*h[0]);
    b[0]=(a[1]-a[0])/h[0]-c[1]*h[0]/3;
    // Calculo de b[k]
    for( int k=0; k<nint-1; k++) b[k+1]=b[k]+(c[k]+c[k+1])*h[k];
    //Impresion de coeficientes
    cout<< " a "<<" b "<<" c "<< " d "<<endl;
    for(int i=0; i<nint; i++) cout<<setprecision(4)<<setw(8)<<a[i]<<" "<<setw(8)<<b[i]<<" "<<setw(8)<<c[i]<<" "<<setw(8)<<d[i]<<endl;
    cout<<endl;
    return 0;
}

template <class T> int
splines_sujetos(Vector<T> &x, Vector<T> &a, Vector<T> &b,
                Vector<T> &c, Vector<T> &d, T& fp0, T& fp1)
{
    cout<<"splines sujetos"<<endl;
    int n=x.size(); //n numero de puntos
    int nint= n-1; //numero de intervalos
    // neq =nint -1; //numero de ecuaciones
    //Creacion vector diagonal superior
    Vector<T> h(nint);
    for( int i=0; i<nint; i++) h[i]=x[i+1]-x[i];
    //Creacion Vector diagonal y de terminos independientes
    Vector<T> u(n), v(n);
    // Valores especiales de u(0), v(0)
    u[0]=2*h[0];
    v[0]= 3*( (a[1]-a[0])/h[0] -fp0);
    for( int i=1; i<nint; i++) u[i]=2*(h[i-1]+h[i]);
    for( int i=0; i<nint-1; i++)
        v[i+1]=3*( (a[i+2]-a[i+1])/h[i+1] - (a[i+1]-a[i])/h[i] );
    //Valores especiales de u(n-1), v(n-1)
    v[nint]=3*(fp1- (a[nint]-a[nint-1])/h[nint-1] );
    u[nint]=2*h[nint-1];
    //Eliminacion de Gauss
    for( int i=0; i<nint; i++)
    {
        u[i+1]=u[i+1]-h[i]*h[i]/u[i];
        v[i+1]=v[i+1]-h[i]*v[i]/u[i];
    }
    //Calculo de c[k] y d[k]
    //Sustitucion hacia atras
    c[nint]=v[nint]/u[nint];
    for( int k=nint-1; k>=0; k--)
    {
        c[k]=(v[k]-h[k]*c[k+1])/u[k];
        d[k]=(c[k+1]-c[k])/(3*h[k]);
    }
    //Calculo de b[k]
    //Valores especiales en x0
    b[0]=fp0;
    for( int k=0; k<nint-1; k++) b[k+1]=b[k]+(c[k]+c[k+1])*h[k];
    //Impresion de coeficientes
    cout<< " a "<<" b "<<" c "<< " d "<<endl;
    for(int i=0; i<nint; i++) cout<<setprecision(4)<<setw(8)<<a[i]<<" "<<
                                    setw(8)<<b[i]<<" "<<setw(8)<<c[i]<<" "<<setw(8)<<d[i]<<endl;
    return 0;
}
template <class T> T
eval_splines(T xx, const Vector<T> &x, const Vector<T> &a,
             const Vector<T> &b, const Vector<T> &c, const Vector<T> &d)
{
    //buscar indice k
    int n=x.size();
    int k=0;
    int kinf=0;
    int ksup=n-1;
    while(ksup-kinf>1)
    {
        k=(ksup+kinf)/2;
        if (x[k]>xx) ksup=k;
        else kinf=k;
    }
    k=kinf;
    // cout<<"xx= "<<xx<<" k= "<<k;
    xx=xx-x[k];
    return a[k]+ xx*( b[k]+ xx*( c[k] +d[k]*xx ) );
}

int main(){
    int nd;
    cout<< "Entrar numero de nodos"<<endl;
    cin>>nd;
    //Construccion de vectores de nodos
    Vector<double> x(nd);
    cout<<"Entrar valores de los nodos"<<endl;
    for (int i=0; i<nd; i++) cin>>x[i];
    double x0=x[0], x1=x[nd-1];
    //Eleccion de la funcion a interpolar
    int ft=-1;
    while((ft!=0)&&(ft!=1)&&(ft!=2)&&(ft!=3))
    {
        cout<<"Entrar funcion: 0- Runge 1- Bernstein 2- sin(x)"
    <<" 3- user "<<endl;
        cin>>ft;
    }
    //Calculo de las ordenadas de interpolacion y de los coeficientes del
    //polinomio
    Vector < double > y(nd);
    for (int i=0; i< nd; i++) y[i]=fun(x[i],ft);
    //Derivadas en los extremos
    double fp0,fp1;
    //Calculo de las derivadas en los extremos
    //Se puede fijar el error de la derivada por extrapolacion
    double hhh=0.00001;
    fp0= (fun(x0+hhh,ft)-fun(x0-hhh,ft))/(2*hhh);
    fp1=(fun(x1+hhh,ft)-fun(x1-hhh,ft))/(2*hhh);
    int nint= nd - 1;
    Vector<double> b(nint), c(nd), d(nint);
    //Impresion de los nodos
    cout<<" x y"<<endl;
    for(int i=0;i<nd;i++)
        cout<<setprecision(4)<<setw(8)<<x[i]<<" "<<setw(8)<<y[i]<<endl;
    cout<<endl<<endl;
    //Splines naturales
    //Calculo de los coeficientes
    Vector < double > a=y;
    splines<double>(x,y,b,c,d);
    //evaluacion del polinomio y funcion en los puntos de dibujo (200)
    int n=200;
    double hh=(x1-x0)/n;
    n++;
    Vector <double> xx(n), yy(n) ,ff(n);
    cout<<" x " << " f(x) "<< " splin(x)" <<endl;
    for (int i = 0; i <n; i++)
    {
        xx[i]=x0+ i*hh;
        yy[i] = eval_splines<double>(xx[i], x, y, b,c,d);
        ff[i] = fun(xx[i],ft);
        if(!(i-(i/5)*5))
            cout<<setw(6)<<xx[i]<<" "<<setprecision(8)<<
                  setw(10)<<ff[i]<<" "<<setw(8)<<yy[i]<<endl;
    }
    cout<<endl;
    //Dibujo del polinomio y la funcion
    Gnuplot g1 = Gnuplot();
    // Seleccionamos escala logaritmica para ciertas funciones
    // if( (ft==0)||(ft==1)) g1.cmd("set logscale y");
    // if( (ft==0)) g1.cmd("set logscale y");
    g1.set_style("lines");
    g1.plot_xy(xx,ff," ");
    g1.set_style("points");
    g1.plot_xy(xx,yy," splines naturales ");
    //Splines sujetos
    //Calculo de los coeficientes
    splines_sujetos<double>(x,y,b,c,d, fp0,fp1);
    cout<<endl<<" x " << " f(x) "<< " splin(x)" <<endl;
    for (int i = 0; i <n; i++)
    {
        xx[i]=x0+ i*hh;
        yy[i] = eval_splines<double>(xx[i], x, y, b,c,d);
        ff[i] = fun(xx[i],ft);
        if(!(i-(i/5)*5))
            cout<<setw(6)<<xx[i]<<" "<<setprecision(8)<<
                  setw(10)<<ff[i]<<" "<<setw(8)<<yy[i]<<endl;
    }
    //Dibujo del polinomio y la funcion
    Gnuplot g2 = Gnuplot();
    // Seleccionamos escala logaritmica para ciertas funciones
    // if( (ft==0)||(ft==1)) g1.cmd("set logscale y");
    // if( (ft==0)) g1.cmd("set logscale y");
    g2.set_style("lines");
    g2.plot_xy(xx,ff," ");
    g2.set_style("points");
    g2.plot_xy(xx,yy,"splines sujetos ");
    //Eleccion del tiempo de permanencia de la grafica en pantalla
    int SLEEP_LGTH=15;
    sleep(SLEEP_LGTH);
    return 0;
}
*/