F*CK beesplines

#include <iostream>
#include <math.h>
#include "stdio.h"

using namespace std;

#define xtype double
#define N 10

xtype _a = 2.477520;
xtype _b = 0.232206;


xtype xs[] = {0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5};
xtype ys[] = {2.78, 3.13, 3.51, 3.94, 4.43, 4.97, 5.58, 6.27, 7.04, 7.91, 8.89};

xtype xsNew[] = {0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 2.00, 2.25, 2.50, 2.75, 3.00, 3.25, 3.50, 3.75, 4.00, 4.25, 4.50, 4.75, 5.00, 5.25, 5.50};
xtype ysNew[] = {2.78, 2.95, 3.13, 3.31, 3.51, 3.72, 3.94, 4.18, 4.43, 4.69, 4.97, 5.27, 5.58, 5.92, 6.27, 6.65, 7.04, 7.47, 7.91, 8.38, 8.89};

xtype f(xtype x) { return _a * exp(_b*x); }

xtype h = xs[1] - xs[0];

typedef struct spline {
    xtype a[N];
    xtype b[N];
    xtype c[N];
    xtype d[N];
} _spline;

typedef struct bispline {
    xtype x[N];
} _bispline;

void solveMatrix(int n, xtype* a, xtype* b, xtype* c, xtype* d, xtype* x) {

    int i;
    xtype alpha[n], beta[n];
    alpha[0] = -c[0] / b[0];
    beta[0] = d[0] / b[0];

    for (i = 1; i < n-1; ++i) {
        alpha[i] = (-c[i]) / (a[i] * alpha[i-1] + b[i]);
        beta[i] = (d[i] - a[i] * beta[i-1]) / (a[i] * alpha[i-1] + b[i]);
    }

    beta[n-1] = (d[n-1]-a[n-1] * beta[n-2]) / (a[n-1] * alpha[n-2] + b[n-1]);
    x[n-1] = beta[n-1];

    for (i = n-2; i >= 0; --i) {
        x[i] = alpha[i] * x[i+1] + beta[i];
    }
}




_spline populateSpline(int n) {

    xtype ak[n-3], bk[n-2], ck[n-3], dk[n-2], ek[n];
    xtype a[n], b[n], c[n], d[n];

    xtype h = xs[1] - xs[0];

    for (int i = 0; i < n - 2; i++) {

        if (i < n-3)
            ak[i] = ck[i] = 1;

        bk[i] = 4;
        dk[i] = 3 / (h * h) * (ys[i+2] - 2*ys[i+1] + ys[i]);
    }

    solveMatrix(n-2, ak, bk, ck, dk, ek);


    c[0] = c[n-1] = 0;
    for (int i = 1; i < n-1; i++)
        c[i] = ek[i-1];

    for (int i = 0; i < n-1; i++) {
        a[i] = ys[i];
        b[i] = (ys[i+1] - ys[i]) / h - (h/3) * (c[i+1] + 2*c[i]);
        d[i] = (c[i+1] - c[i]) / (3*h);
    }


    _spline spZ;

    for (int i = 0; i < n; i++) {
        spZ.a[i] = a[i];
        spZ.b[i] = b[i];
        spZ.c[i] = c[i];
        spZ.d[i] = d[i];
    }
    return spZ;
}

xtype countSplineAt(_spline spline, xtype x) {

    //S_j[x] = a_j + b_j(x-xj) + c_j(x-xj)^2 + d_j(x-xj)^3
    int i;
    for (i = 0; i < N; ++i)
        if (x >= xs[i] && x <= xs[i+1])
            break;

    if (N == i) i--;

    //printf("cSA: %.6f\n", xsNew[i]);
    //printf("[%d] %.6f %.6f \n",i, x, xsNew[i]);

    return  spline.a[i] +
            spline.b[i]*(x - xs[i]) +
            spline.c[i]*(x - xs[i])*(x - xs[i]) +
            spline.d[i]*(x - xs[i])*(x - xs[i])*(x - xs[i]);
}


xtype nFunc(int q, int k, xtype t) {

    if (1 == q) {
        if (t >= xs[k] && t <= xs[k+1])
            return 1;
        else
            return 0;
    }

    xtype n1 = nFunc(q-1, k, t) * (t - xs[k]) / (xs[k+q-1] - xs[k]);
    xtype n2 = nFunc(q-1, k+1, t) * (xs[k+q] - t) / (xs[k+q] - xs[k+1]);

    return n1 + n2;
}



bispline populateBiSpline(int n) {

    xtype ak[n], bk[n], ck[n], dk[n], ek[n+2];
    xtype a[n], b[n], c[n], d[n];

    for (int i = 1; i < n - 2; i++) {
        a[i] = 1;
        b[i] = 4;
        c[i] = 1;
        d[i] = 6*ys[i];
    }

    a[0] = 0;
    b[0] = 1;
    c[0] = 0;
    d[0] = ys[0];

    a[n-3] = 0;
    b[n-3] = 1;
    c[n-3] = 0;
    d[n-3] = ys[N-3];

    for (int i = 0; i < n; i++) {
        d[i] /= h*h*h;
    }

    solveMatrix(n-2, a, b, c, d, ek);


    _bispline bsp;

    for (int i = 1; i < n-1; i++)
        bsp.x[i] = ek[i-1];

    bsp.x[0] = 2*xs[1] - xs[2];
    bsp.x[n-2] = 2*xs[n-2] - xs[n-3];

    return bsp;
}

xtype baseX(int k, xtype x, xtype kh) {
    if (0 == k) {

        if (xs[0] + 2*h <= x) return 0;
        if (xs[0] + h <= x && x < xs[0] + 2*h) return pow(2*h + xs[0] - kh - (x - kh), 3)/6;
        if (xs[0] <= x && x < xs[0] + h)
            return 2*pow(h, 3)/3 - pow(-xs[0] + kh + (x - kh), 2)*(2*h + xs[0] - kh - (x - kh))/2;
        if (xs[0] - h <= x && x < xs[0])
            return 2*pow(h, 3)/3 - pow(-xs[0] + kh + (x - kh), 2)*(2*h - xs[0] + kh + (x - kh))/2;
        if (xs[0] - 2*h <= x && x < xs[0] - h) return pow(2*h - xs[0] + kh + (x - kh), 3)/6;
        if (x < xs[0] - 2*h) return 0;

    } else {
        return baseX(0, x - k*h, -k*h);
    }

}

xtype countBiSplineAt(_bispline bsp, xtype val) {

    int i;
    for (i = 0; i < N; ++i)
        if (xs[i] <= val && val <= xs[i+1])
            break;

    if (N == i) i--;

    //printf("%.2f<=%.2f<=%.2f\n", xs[i], val, xs[i+1]);
    if (xs[i] <= val && val <= xs[i+1]) {
        return  bsp.x[i + 0] * baseX(i - 1, val, 0) +
                bsp.x[i + 1] * baseX(i + 0, val, 0) +
                bsp.x[i + 2] * baseX(i + 1, val, 0) +
                bsp.x[i + 3] * baseX(i + 2, val, 0);
    } else {
        return -1;
    }
}


int main()
{
    printf("populating spline\n");
    _spline sp = populateSpline(N+1);
    _bispline bs = populateBiSpline(N+1);


    for ( int i = 0; i <= 2*N; i++) {

        xtype sVal = countSplineAt(sp, xsNew[i]);

        xtype bsVal = countBiSplineAt(bs, xsNew[i]);

        printf("x[%d]=%.2f y=%.5f s=%.5f d=%.5f b=%.5f d=%.5f\n",
               i, xsNew[i], ysNew[i], sVal, fabs(ysNew[i] - sVal), bsVal, fabs(ysNew[i] - bsVal));

        //printf("x[%d]=%.2f y[%d]=%.6f bs[%d]=%.6f\n",
        //        i, xsNew[i], i, ysNew[i], i, bsVal);
    }

/*  xtype s = 0;
    int k = 3;
    int q = 4;


    for ( int k = 0; k < N; k++) {

        s += nFunc(q, k, 0.1) * ys[k];
    }
    printf("s=%.6f", s);*/

    return 0;
}