num meth splines

#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+3];
} _bispline;

typedef struct akame {
    xtype a[N+4];
    xtype b[N+4];
    xtype c[N+4];
    xtype d[N+4];
} _akame;

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

void solveMatrixNew(INT len, xtype* a, xtype* b, xtype* c, xtype* d, xtype* x) {

    INT SIZE = len;
    INT n = SIZE - 1;
    INT i;

    xtype alpha[SIZE], beta[SIZE];

    alpha[0] = -c[0] / b[0];
    beta[0] = d[0] / b[0];

    FOR (i = 1; i <= n; 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];

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




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 len) {

    INT SIZE = len;
    INT n = SIZE - 1;

    xtype a[SIZE], b[SIZE], c[SIZE], d[SIZE], e[SIZE+2];

    FOR (INT i = 1; i < n; 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] = 0;
    b[n] = 1;
    c[n] = 0;
    d[n] = ys[n-1];

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

    solveMatrixNew(len, a, b, c, d, e);


    _bispline bsp;


    FOR (INT i = 0; i < SIZE+2; i++)
        bsp.x[i] = e[i];

    RETURN bsp;
}

xtype baseX(INT k, xtype x, xtype kh) {
    IF (0 == k) {

        xtype xz = xs[0];

        IF (xz + 2*h <= x) RETURN 0;
        IF (xz + h <= x && x < xz + 2*h) RETURN pow(2*h + xz - kh - (x - kh), 3)/6;
        IF (xz <= x && x < xz + h)
            RETURN 2*pow(h, 3)/3 - pow(-xz + kh + (x - kh), 2)*(2*h + xz - kh - (x - kh))/2;
        IF (xz - h <= x && x < xz)
            RETURN 2*pow(h, 3)/3 - pow(-xz + kh + (x - kh), 2)*(2*h - xz + kh + (x - kh))/2;
        IF (xz - 2*h <= x && x < xz - h) RETURN pow(2*h - xz + kh + (x - kh), 3)/6;
        IF (x < xz - 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+2; ++i)
        IF (xs[i] <= val && val <= xs[i+1])
            break;

    IF (N+2 == 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;
    //}
}


akame populateAkame(INT n) {

    INT i;
    xtype m[n + 4];
    //-2;-1;    0.. n-1    ;n;n+1
    FOR (i = 2; i < n + 2; i++) {
        m[i] = (ys[i+1-2] - ys[i-2]) / (xs[i+1-2] - xs[i-2]);
    }

    m[0] = 3*m[2] - 2*m[3]; //-2=0-1
    m[1] = 2*m[2] - 2*m[3];
    m[n+2] = 2*m[n+1] - 2*m[n];
    m[n+3] = 3*m[n];

    xtype ne, alpha_i, h_i;
    xtype t_l[n], t_r[n];

    FOR (i = 2; i < n + 2; i++) {
        ne = fabs(m[i+1] - m[i]) + fabs(m[i-1] - m[i-2]);

        IF (ne > 0) {
            alpha_i = fabs(m[i-1] - m[i-2]) / ne;
            t_l[i] = m[i-1] + alpha_i*(m[i] - m[i-1]);
            t_r[i] = t_l[i];
        } ELSE {
            t_l[i] = m[i-1];
            t_r[i] = m[i];
        }
    }

    akame ga_spline;
    FOR (i = 2; i < n+2; i++) {
        h_i = xs[i+1-2] - xs[i-2];

        ga_spline.a[i] = ys[i-2];
        ga_spline.b[i] = t_r[i];
        ga_spline.c[i] = (3*m[i] - 2*t_r[i] - t_l[i+1]) / h_i;
        ga_spline.d[i] = (t_r[i] + t_l[i+1] - 2*m[i]) / (h_i*h_i);
    }

    RETURN ga_spline;
}

xtype countAkameAt(akame ak, xtype x, INT i) {

    FOR (i = 0; i < N; ++i)
        IF (x >= xs[i] && x <= xs[i+1])
            break;

    RETURN  ak.a[i+2] +
            ak.b[i+2]*(x - xs[i]) +
            ak.c[i+2]*(x - xs[i])*(x - xs[i]) +
            ak.d[i+2]*(x - xs[i])*(x - xs[i])*(x - xs[i]);
}

INT f1() {
    CHAR ch, ch1 = 'a', ch2 = 'b', ch3 = 'c';
    ch = ch1 + ch2 + ch3;
    //f1 = INT(ch);
    //RETURN f1;
}

INT SUM(INT START, INT END) {
    IF (END == START)
        RETURN 1;
    ELSE
        RETURN END + SUM(START, END - 1);
}

INT main()
{
    printf("%d\n", SUM(1,3));

    printf("populating Spline\n");
    spline sp = populateSpline(N+1);

    printf("populating Bispline\n");
    bispline bs = populateBiSpline(N+2);

    printf("populating Akame ga Spline\n");
    akame AS = populateAkame(N);

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

        xtype sVal = countSplineAt(sp, xsNew[i]);
        xtype bsVal = countBiSplineAt(bs, xsNew[i]);
        xtype asVal = countAkameAt(AS, xsNew[i], i);

        printf("x[%d]=%.2f y=%.5f \ts=%.5f d=%.5f \tbi=%.5f d=%.5f \ta=%.5f d=%.5f\n",
               i, xsNew[i], ysNew[i],
               sVal, fabs(ysNew[i] - sVal),
               bsVal, fabs(ysNew[i] - bsVal),
               asVal, fabs(ysNew[i] - asVal)
               );
    }
    RETURN 0;
}