// Given a bunch of evenly spaced points, interpolate at tenths between the
// middle two points, using 1st, 3rd, 5th, 7th, 9th degree Lagrange
// polynomials.  Pretend given points are on odd x, with 0 in the center.

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

#include <stdio.h>

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

void third(double a, double b, double c, double d, double* coef)
{
    double p1 = c+b, m1 = c-b;
    double p2 = d+a, m2 = d-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 fifth(
    double a,
    double b,
    double c,
    double d,
    double e,
    double f,
    double *coef)
{
    double p1 = d+c, m1 = d-c;
    double p2 = e+b, m2 = e-b;
    double p3 = f+a, m3 = f-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 seventh(
    double a,
    double b,
    double c,
    double d,
    double e,
    double f,
    double g,
    double h,
    double* coef)
{
    double p1 = e+d, m1 = e-d;
    double p2 = f+c, m2 = f-c;
    double p3 = g+b, m3 = g-b;
    double p4 = h+a, m4 = h-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 ninth(
    double a,
    double b,
    double c,
    double d,
    double e,
    double f,
    double g,
    double h,
    double i,
    double j,
    double* coef)
{
    double p1 = f+e, m1 = f-e;
    double p2 = g+d, m2 = g-d;
    double p3 = h+c, m3 = h-c;
    double p4 = i+b, m4 = i-b;
    double p5 = j+a, m5 = j-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;
}


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];
    first(a[4],a[5],&coef[1][0]);
    third(a[3],a[4],a[5],a[6],&coef[3][0]);
    fifth(a[2],a[3],a[4],a[5],a[6],a[7],&coef[5][0]);
    seventh(a[1],a[2],a[3],a[4],a[5],a[6],a[7],a[8],&coef[7][0]);
    ninth(a[0],a[1],a[2],a[3],a[4],a[5],a[6],a[7],a[8],a[9],&coef[9][0]);

    // print the coefficients
    for (int s=1; s<10; s+=2)
    {
        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=1; s<10; s+=2)
        {
            double r = 0;
            double xn = 1.0;
            for (int t=0; t<=s; ++t)
            {
                r += coef[s][t]*xn;
                xn *= x;
            }
            printf("%d: %.20f  ", s, r);
        }
        printf("\n");
    }
    return 0;
}