// ispline.cpp 
#include <math.h>

typedef double flt;
typedef flt flt4[4];
typedef flt flt7[7];

inline flt sqr(flt x) { return x*x; }

bool ispline3_build( // ispline3_build - строит аппроксимирующий сплайн 3-го порядка
    int n,   // кол-во точек
    const flt *x,const flt *f,const flt *df, // [n]
    flt sm,  // критерий гладкости 
    flt4 *c, // [n]   результат (коэф. для сплайна)
    flt7 *t  // [n+2] временные данные
) 
// n - кол-во точек
// x,f,df - x-сетка f зн ф-и, df погр f (n штук)
// c - результат (коэф. для сплайна)
// t - рабочий массив (n+2 элементов)                         
// sm - критерий гладкости (0..inf)
{
    flt g=0;
    t[0  ][0]=0; t[1  ][0]=0;
    t[n  ][2]=0; t[n+1][2]=0;
    t[n  ][1]=0;
    t[0  ][5]=0; t[n  ][5]=0;
    t[1  ][5]=0; t[n+1][5]=0;
    flt p=0,h=x[1]-x[0]; if (h<=0) return false;
    flt f2=-sm, ff=(f[1]-f[0])/h;
    {for(int i=2;i<n;i++) {
        g=h; h=x[i]-x[i-1]; if (h<=0) return false;
        flt ldh=1/h;
        flt e=ff; ff=(f[i]-f[i-1])*ldh;
        c[i][0]=ff-e;
        t[i][3]=(g+h)*2.0/3.0;
        t[i][4]=h/3.0;
        t[i][2]=df[i-2]/g;
        t[i][0]=df[i]*ldh;
        t[i][1]=-df[i-1]*(1.0/g+ldh);
    }}
    {for(int i=2;i<n;i++) {
        c[i-1][1]=sqr(t[i][0])+sqr(t[i][1])+sqr(t[i][2]);
        c[i-1][2]=t[i][0]*t[i+1][1]+t[i][1]*t[i+1][2];
        c[i-1][3]=t[i][0]*t[i+2][2];
    }}
    for(;;) {
        {for(int i=2;i<n;i++) {
          t[i-1][1]=ff*t[i-1][0];
          t[i-2][2]= g*t[i-2][0];
          t[i  ][0]=1.0/(p*c[i-1][1]+t[i][3]-ff*t[i-1][1]-g*t[i-2][2]);
          t[i  ][5]=c[i][0]-t[i-1][1]*t[i-1][5]-t[i-2][2]*t[i-2][5];
          ff=p*c[i-1][2]+t[i][4]-h*t[i-1][1];
          g=h; h=c[i-1][3]*p;
        }}
        {for(int i=2;i<n;i++) {
          int j=n+1-i;
          t[j][5]=t[j][0]*t[j][5]-t[j][1]*t[j+1][5]-t[j][2]*t[j+2][5];
        }}
        flt e=0; h=0;
        {for(int i=1;i<n;i++) {
          g=h; h=(t[i+1][5]-t[i][5])/(x[i]-x[i-1]);
          flt hmg=h-g;
          t[i][6]=hmg*sqr(df[i-1]);
          e=e+t[i][6]*hmg;
        }}
        g=-h*sqr(df[n-1]);
        t[n][6]=g;
        e=e-g*h;
        g=f2;
        f2=e*p*p;
        if (f2>=sm || f2<=g)  break;
        ff=0;
        h=(t[2][6]-t[1][6])/(x[1]-x[0]);
        {for(int i=2;i<n;i++) {
          g=h; h=(t[i+1][6]-t[i][6])/(x[i]-x[i-1]);
          g=h-g-t[i-1][1]*t[i-1][0]-t[i-2][2]*t[i-2][0];
          ff=ff+g*t[i][0]*g;
          t[i][0]=g;
        }}
        h=e-p*ff;
        if (h<=0) break;
        p=p+(sm-f2)/((sqrt(sm/e)+p)*h);
    }
    {for(int i=0;i<n-1;i++) {
        c[i][0]=f[i]-p*t[i+1][6];
        c[i][2]=t[i+1][5];
        t[i][0]=c[i][0];
    }}
    t[n-1][0]=f[n-1]-p*t[n][6];
    c[n-1][0]=t[n-1][0];
    {for(int i=1;i<n;i++) {
        h=x[i]-x[i-1];
        c[i-1][3]=(t[i+1][5]-c[i-1][2])/h/3.0;
        c[i-1][1]=(t[i  ][0]-c[i-1][0])/h - (h*c[i-1][3]+c[i-1][2])*h;
    }}
    return true;
}

inline flt spln(flt t,const flt4 &s) { return ((s[3]*t+s[2])*t+s[1])*t+s[0]; } // poly3

flt spline(flt x,int n,const flt *xn,const flt4* c) 
// x - расчетная точка
// xn - сетка : xn[i]<xn[i+1]
// с - коэф. сплайна (строятся с помощью ispline3_build)
{
    int i=0,j=n-1;
    if (x<xn[i  ]) j=i+1;
    if (x>xn[j-1]) i=j-1;
    while ((j-i)>1) {
        int m=(i+j)>>1;
        if (x<xn[m]) j=m; else i=m;
    }
    return spln(x-xn[i],c[i]);
}             

/*
    void example() {
        enum { n=4 };
        flt  x[n]={1,2,3,4};
        flt  y[n]={4,2,0,4};
        flt dy[n]={0.1,1,3,0.1};
        flt7 tmp[n+2]; flt4 c[n];
        flt sm=10.0;
        ispline3_build(n,x,y,dy,sm,c,tmp);
        flt xo=2.5;
        flt fo=spline(xo,n,x,c);
    }
*/