cubic spline accepted

#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); }

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

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]);
}

int main()
{
    int i;

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

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

    for (i = 1; i < 2*N; 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+1);



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

        //printf("in: %.6f\n", xsNew[i]);
        xtype sVal = countSplineAt(sp, xsNew[i]);

        //xtype sValMid = countSplineAt(sp, xsNew[i] + 0.125);

        printf("x[%d]=%.6f y[%d]=%.6f s[%d]=%.6f d=%.6f\n",
               i, xsNew[i], i, ysNew[i], i, sVal, fabs(ysNew[i] - sVal));
        //printf("x[%d]=%.6f                s[%d]=%.6f\n",
        //        i, xsNew[i] + 0.125, i, sValMid);

        //printf ("(%.3f; %.3f) ", xsNew[i], sVal);
        //printf ("(%.3f; %.3f) ", xsNew[i] + 0.125, sValMid);
    }

    return 0;
}



/*
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;
}*/