chmet3 cubic beespline

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

using namespace std;

#define xtype float
#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};
xtype ys[] = {3.02, 3.21, 2.95, 4.06, 4.03, 5.39, 5.97, 6.51, 6.77, 7.79};

xtype xsNew[N*2];
xtype ysNew[N*2];

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

typedef struct spline {
    xtype a[N*2-3];
    xtype b[N*2-2];
    xtype c[N*2-3];
    xtype d[N*2-2];
} _spline;

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

// |b a . .|  |d|  |x|
// |c b a .|  |d|  |x|
// |. c b a|  |d|  |x|
// |. . c b|  |d|  |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 a[n-3], b[n-2], c[n-3], d[n-2], e[n];

    xtype h = (xsNew[n-1] - xsNew[0]) / n;

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

        if (i < n-3)
            a[i] = c[i] = 1;

        b[i] = 4;
        d[i] = 3 / (h * h) * (ysNew[i+2] - 2*ysNew[i+1] + ysNew[i]);
    }

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

    c[0] = c[n-1] = 0;

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

    for (int i = 0; i < n-1; ++i) {
        a[i] = ysNew[i];
        b[i] = (ysNew[i+1] - ysNew[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];
        printf("[%d] a:%.3f b:%.3f c:%.3f d:%.3f\n",
                i, a[i], b[i], c[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-1; ++i)
        if (x >= xsNew[i] && x <= xsNew[i+1])
            break;

    if (N-1 == i) i--;

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

xtype nFunc(int n, int i, xtype x) {

    if (x >= xsNew[i] && x <= xsNew[i+1])
        return 1;
    else
        return 0;

    xtype n1 = nFunc(n-1, i, x) * (x - xsNew[i]) / (xsNew[i+n] - xsNew[i]);
    xtype n2 = nFunc(n-1, i+1, x) * (xsNew[i+n+1] - x) / (xsNew[i+n+1] - xsNew[i+1]);

    return n1 + n2;
}


xtype dFunc(int k, int i, xtype x) {

    xtype alpha = (x - xsNew[i]) / (xsNew[i+N*2+1-k] - xsNew[i]);

    xtype d1 = dFunc(k-1, i-1, x);
    xtype d2 = dFunc(k-1, i, x);

    return (1 - alpha)*d1 + alpha*d2;
}

int main()
{
    int i;

    xsNew[0] = 0.25;
    ysNew[0] = f(xsNew[0]);

    printf("increasing frequency of X-values\n");
    printf("[0](%.2f, %f)\n", xsNew[0], ysNew[0]);

    for (i = 1; i < 20; i++) {

        xsNew[i] = xsNew[i-1] + 0.25;
        ysNew[i] = f(xsNew[i]);

        printf("[%d](%.2f, %f)\n", i, xsNew[i], ysNew[i]);
    }

    printf("populating spline\n");
    _spline sp = populateSpline(N*2);


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

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

        printf("y[%d]=%f s[%d]=%.16f d=%f\n",
                i, ysNew[i], i, sVal, fabs(ysNew[i] - sVal));
    }

    return 0;
}