// Given a bunch of evenly spaced points, interpolate between the
// middle two points, as if the points are spaced 2 apart and 0 were
// in the middle of the range.
// Input: values points 0..n
// Input: coefficients 0..n of the interpolation polynomial

// typedef long double ld;  sscanf("%llf"); printf("%Lf");

#include <stdio.h>

void zeroth(double* coef)
{
    // coef[0] is already equal to coef[0]
    return;
}

void first(double* coef)
{
    double a,b;
    a = coef[0], b = coef[1];
    double p1 = b+a, m1 = b-a;
    coef[0] = p1 / 2.0;
    coef[1] = m1 / 2.0;
}

void second(double* coef)
{
    double c = coef[1];
    double a,b;
    a = coef[0], b = coef[2];
    double p1 = b+a, m1 = b-a;
    coef[0] = c;
    coef[1] = m1 / 4;
    coef[2] = (-2*c + p1) / 8;
}

void third(double* coef)
{
    double a,b;
    a = coef[1], b = coef[2];
    double p1 = b+a, m1 = b-a;
    a = coef[0], b = coef[3];
    double p2 = b+a, m2 = b-a;
    coef[0] = (9*p1 - p2) / 16;
    coef[1] = (27*m1 - m2) / 48;
    coef[2] = (-p1 + p2) / 16;
    coef[3] = (-3*m1 + m2) / 48;
}

void fourth(double* coef)
{
    double c = coef[2];
    double a,b;
    a = coef[1], b = coef[3];
    double p1 = b+a, m1 = b-a;
    a = coef[0], b = coef[4];
    double p2 = b+a, m2 = b-a;
    coef[0] = c;
    coef[1] = (8*m1 - m2) / 24;
    coef[2] = (-30*c + 16*p1 - p2) / 96;
    coef[3] = (-2*m1 + m2) / 96;
    coef[4] = (6*c - 4*p1 + p2) / 384;
}

void fifth(double *coef)
{
    double a,b;
    a = coef[2], b = coef[3];
    double p1 = b+a, m1 = b-a;
    a = coef[1], b = coef[4];
    double p2 = b+a, m2 = b-a;
    a = coef[0], b = coef[5];
    double p3 = b+a, m3 = b-a;
    coef[0] = (150*p1 - 25*p2 + 3*p3) / 256;
    coef[1] = (2250*m1 - 125*m2 + 9*m3) / 3840;
    coef[2] = (-34*p1 + 39*p2 - 5*p3) / 384;
    coef[3] = (-34*m1 + 13*m2 - m3) / 384;
    coef[4] = (2*p1 - 3*p2 + p3) / 768;
    coef[5] = (10*m1 - 5*m2 + m3) / 3840;
}

void sixth(double* coef)
{
    double c = coef[3];
    double a,b;
    a = coef[2], b = coef[4];
    double p1 = b+a, m1 = b-a;
    a = coef[1], b = coef[5];
    double p2 = b+a, m2 = b-a;
    a = coef[0], b = coef[6];
    double p3 = b+a, m3 = b-a;
    coef[0] = c;
    coef[1] = (45*m1 - 9*m2 + m3) / 120;
    coef[2] = (-490*c + 270*p1 - 27*p2 + 2*p3) / 1440;
    coef[3] = (-13*m1 + 8*m2 - m3) / 384;
    coef[4] = (56*c - 39*p1 + 12*p2 - p3) / 2304;
    coef[5] = (5*m1 - 4*m2 + m3) / 7680;
    coef[6] = (-20*c + 15*p1 - 6*p2 + p3) / 46080;
}

void seventh(double* coef)
{
    double a,b;
    a = coef[3], b = coef[4];
    double p1 = b+a, m1 = b-a;
    a = coef[2], b = coef[5];
    double p2 = b+a, m2 = b-a;
    a = coef[1], b = coef[6];
    double p3 = b+a, m3 = b-a;
    a = coef[0], b = coef[7];
    double p4 = b+a, m4 = b-a;
    coef[0] = (385875*p1 - 77175*p2 + 15435*p3 - 1575*p4) / 645120;
    coef[1] = (385875*m1 - 25725*m2 + 3087*m3 - 225*m4) / 645120;
    coef[2] = (-66185*p1 + 81837*p2 - 17465*p3 + 1813*p4) / 645120;
    coef[3] = (-66185*m1 + 27279*m2 - 3493*m3 + 259*m4) / 645120;
    coef[4] = (2905*p1 - 4725*p2 + 2065*p3 - 245*p4) / 645120;
    coef[5] = (2905*m1 - 1575*m2 + 413*m3 - 35*m4) / 645120;
    coef[6] = (-5*p1 + 9*p2 - 5*p3 + p4) / 92160;
    coef[7] = (-35*m1 + 21*m2 - 7*m3 + m4) / 645120;
}

void eighth(double* coef)
{
    double c = coef[4];
    double a,b;
    a = coef[3], b = coef[5];
    double p1 = b+a, m1 = b-a;
    a = coef[2], b = coef[6];
    double p2 = b+a, m2 = b-a;
    a = coef[1], b = coef[7];
    double p3 = b+a, m3 = b-a;
    a = coef[0], b = coef[8];
    double p4 = b+a, m4 = b-a;
    coef[0] = c;
    coef[1] = (672*m1 - 168*m2 + 32*m3 - 3*m4) / 1680;
    coef[2] = (-14350*c + 8064*p1 - 1008*p2 + 128*p3 - 9*p4) / 40320;
    coef[3] = (-488*m1 + 338*m2 - 72*m3 + 7*m4) / 11520;
    coef[4] = (2730*c - 1952*p1 + 676*p2 - 96*p3 + 7*p4) / 92160;
    coef[5] = (29*m1 - 26*m2 + 9*m3 - m4) / 23040;
    coef[6] = (-150*c + 116*p1 - 52*p2 + 12*p3 - p4) / 184320;
    coef[7] = (-14*m1 + 14*m2 - 6*m3 + m4) / 1290240;
    coef[8] = (70*c - 56*p1 + 28*p2 - 8*p3 + p4) / 10321920;
}

void ninthy(double* coef)
{
    double a,b;
    a = coef[4], b = coef[5];
    double p1 = b+a, m1 = b-a;
    a = coef[3], b = coef[6];
    double p2 = b+a, m2 = b-a;
    a = coef[2], b = coef[7];
    double p3 = b+a, m3 = b-a;
    a = coef[1], b = coef[8];
    double p4 = b+a, m4 = b-a;
    a = coef[0], b = coef[9];
    double p5 = b+a, m5 = b-a;
    coef[0] = (99225*p5 - 1148175*p4 + 6429780*p3 - 25004700*p2 + 112521150*p1) / 185794560;
    coef[1] = (11025*m5 - 164025*m4 + 1285956*m3 - 8334900*m2 + 112521150*m1) / 185794560;
    coef[2] = (-116244*p5 + 1335852*p4 - 7354800*p3 + 26823888*p2 - 20688696*p1) / 185794560;
    coef[3] = (-12916*m5 + 190836*m4 - 1470960*m3 + 8941296*m2 - 20688696*m1) / 185794560;
    coef[4] = (17766*p5 - 194922*p4 + 950040*p3 - 1858248*p2 + 1085364*p1) / 185794560;
    coef[5] = (1974*m5 - 27846*m4 + 190008*m3 - 619416*m2 + 1085364*m1) / 185794560;
    coef[6] = (-756*p5 + 7308*p4 - 25200*p3 + 39312*p2 - 20664*p1) / 185794560;
    coef[7] = (-84*m5 + 1044*m4 - 5040*m3 + 13104*m2 - 20664*m1) / 185794560;
    coef[8] = (9*p5 - 63*p4 + 180*p3 - 252*p2 + 126*p1) / 185794560;
    coef[9] = (m5 - 9*m4 + 36*m3 - 84*m2 + 126*m1) / 185794560;
}


void ninth(double* coef)
{
    double a,b;
    a = coef[4], b = coef[5];
    double p1 = b+a, m1 = b-a;
    a = coef[3], b = coef[6];
    double p2 = b+a, m2 = b-a;
    a = coef[2], b = coef[7];
    double p3 = b+a, m3 = b-a;
    a = coef[1], b = coef[8];
    double p4 = b+a, m4 = b-a;
    a = coef[0], b = coef[9];
    double p5 = b+a, m5 = b-a;
    coef[0] = (39690*p1 - 8820*p2 + 2268*p3 - 405*p4 + 35*p5) / 65536;
    coef[1] = (11025*m5 - 164025*m4 + 1285956*m3 - 8334900*m2 + 112521150*m1) / 185794560;
    coef[2] = (-116244*p5 + 1335852*p4 - 7354800*p3 + 26823888*p2 - 20688696*p1) / 185794560;
    coef[3] = (-12916*m5 + 190836*m4 - 1470960*m3 + 8941296*m2 - 20688696*m1) / 185794560;
    coef[4] = (17766*p5 - 194922*p4 + 950040*p3 - 1858248*p2 + 1085364*p1) / 185794560;
    coef[5] = (1974*m5 - 27846*m4 + 190008*m3 - 619416*m2 + 1085364*m1) / 185794560;
    coef[6] = (-756*p5 + 7308*p4 - 25200*p3 + 39312*p2 - 20664*p1) / 185794560;
    coef[7] = (-84*m5 + 1044*m4 - 5040*m3 + 13104*m2 - 20664*m1) / 185794560;
    coef[8] = (9*p5 - 63*p4 + 180*p3 - 252*p2 + 126*p1) / 185794560;
    coef[9] = (m5 - 9*m4 + 36*m3 - 84*m2 + 126*m1) / 185794560;
}


int main(int argc, char**argv)
{
    // get 10 consecutive values
    double a[10];
    if (argc != 11)
    {
        printf("give ten doubles as arguments");
        return 1;
    }
    for (int i=0; i<10; ++i)
    {
        sscanf(argv[i+1], "%lf", &a[i]);
        printf("%d:%f ");
    }
    printf("\n");

    // find the coefficients for Lagrange polynomials of odd degrees
    double coef[10][10];
    for (int i=0; i<=1; ++i) coef[1][i] = a[i+4];
    first(&coef[1][0]);
    for (int i=0; i<=3; ++i) coef[3][i] = a[i+3];
    third(&coef[3][0]);
    for (int i=0; i<=5; ++i) coef[5][i] = a[i+2];
    fifth(&coef[5][0]);
    for (int i=0; i<=7; ++i) coef[7][i] = a[i+1];
    seventh(&coef[7][0]);
    for (int i=0; i<=9; ++i) coef[9][i] = a[i+0];
    ninth(&coef[9][0]);

    // find the coefficients for Lagrange polynomials of even degrees
    for (int i=0; i<=0; ++i) coef[0][i] = a[i+4];
    zeroth(&coef[0][0]);
    for (int i=0; i<=2; ++i) coef[2][i] = a[i+3];
    second(&coef[2][0]);
    for (int i=0; i<=4; ++i) coef[4][i] = a[i+2];
    fourth(&coef[4][0]);
    for (int i=0; i<=6; ++i) coef[6][i] = a[i+1];
    sixth(&coef[6][0]);
    for (int i=0; i<=8; ++i) coef[8][i] = a[i+0];
    eighth(&coef[8][0]);
    
    // print the coefficients
    for (int s=0; s<10; s++)
    {
        printf("%d: ");
        for (int t=0; t<=s; ++t)
            printf(" %g", coef[s][t]);
        printf("\n");
    }

    // do some interpolations with the polynomials
    for (double x=-1.0; x<1.00001; x+=0.5)
    {
        printf("x=%f : ", x);
        for (int s=0; s<10; s++)
        {
            double r = 0;
            double yn = 1.0;
            double y = (s&1) ? x : x+1;
            for (int t=0; t<=s; ++t)
            {
                r += coef[s][t]*yn;
                yn *= y;
            }
            printf("%d: %.16f  ", s, r);
        }
        printf("\n");
    }
    return 0;
}