cubic spline eng wiki algo

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

using namespace std;

#define xtype float
#define n 20

xtype _a = 2.477520;
xtype _b = 0.232206;

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


typedef struct spline {
    xtype a[n + 1];
    xtype b[n];
    xtype c[n + 1];
    xtype d[n];
    xtype x[n];
} _spline;

_spline populateSpline(xtype x[n], xtype y[n]) {
    int i;

    xtype a[n+1], b[n], d[n];
    xtype h[n], alpha[n];
    xtype c[n+1], l[n+1], mu[n+1], z[n+1];


    for (i = 0; i <= n; i++)
        a[i] = y[i];

    for (i = 0; i < n; i++)
        h[i] = x[i+1] - x[i];


    for (i = 1; i < n; i++)
        alpha[i] = (a[i+1] - a[i]) * 3/h[i] -
                        (a[i] - a[i-1]) * 3/h[i-1];


    l[0] = 1, mu[0] = z[0] = 0;

    for (i = 1; i < n-1; i++) {

        l[i] = 2*(x[i+1] - x[i-1]) - h[i-1]*mu[i-1];
        mu[i] = h[i] / l[i];
        z[i] = (alpha[i] - h[i-1]*z[i-1]) / l[i];

        printf("l:%.3f mu:%.3f z:%.3f\n",
                l[i], mu[i], z[i]);
    }



    l[(int)n] = 1;
    z[(int)n] = c[(int)n] = 0;

    for (i = n - 1; i >= 0; --i) {
        c[i] = z[i] - mu[i]*c[i+1];

        b[i] = (a[i+1] - a[i]) / h[i] -
                    h[i]*(c[i+1] + 2*c[i]) / 3;

        d[i] = (c[i+1] - c[i]) / (3 * h[i]);
        printf("b:%.3f c:%.3f d:%.3f\n",
                b[i], c[i], d[i]);
    }


    _spline spZ;
    for (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];
        spZ.x[i] = x[i];
        printf("a:%.3f b:%.3f c:%.3f d:%.3f x:%.3f\n",
                a[i], b[i], c[i], d[i], x[i]);
    }

    return spZ;
}

xtype countSplineAt(_spline spline, int i, xtype x) {

    //S_j[x] = a_j + b_j(x-xj) + c_j(x-xj)^2 + d_j(x-xj)^3

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




int main() {
    int i;
    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];
    xtype ysNew[n];
    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]);
    }
    int iAmount = i;

    printf("populating spline\n");
    _spline spline = populateSpline(xsNew, ysNew);


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

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

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






    return 0;
}