#include <string.h>
#include "biglagrange.h"

// implement the function to interpolate here
BigFloat F(const BigFloat& x, void* context)
{
    BigFloat r(x);
    return r.Sin();
}

int main(int argc, char** argv)
{
    // interpret the arguments
    if (argc != 8 && argc != 11)
    {
        fprintf(stderr, "usage: interpolate name sn sd fn fd nPieces nTerms [--integeral cn cd]\n");
        fprintf(stderr, "start sn/sd to fn/fd, in nPieces equal-sized ranges, with nTerms terms per range\n");
        fprintf(stderr, "if --integral is used, F(start) = cn/cd\n");
        return 1;
    }
    const char* name = argv[1];
    s8 sn, sd, fn, fd, nPieces, nTerms, cn, cd;
    sscanf(argv[2], "%lld", &sn);
    sscanf(argv[3], "%lld", &sd);
    BigFloat start(sn);
    start.Div(sd);
    
    sscanf(argv[4], "%lld", &fn);
    sscanf(argv[5], "%lld", &fd);
    BigFloat finish(fn);
    finish.Div(fd);
    
    sscanf(argv[6], "%lld", &nPieces);
    if (nPieces <= 0)
    {
        fprintf(stderr, "nPieces (number of pieces in the range) must be positive\n");
        return 1;
    }

    BigFloat range(finish);
    range.Sub(start).Div(nPieces);

    sscanf(argv[7], "%lld", &nTerms);
    if (nTerms <= 0)
    {
        fprintf(stderr, "nTerms (number of terms per polynomial) must be positive\n");
        return 1;
    }

    bool integral = false;
    if (argc == 11)
    {
        if (strcmp(argv[8], "--integral"))
        {
            fprintf(stderr, "last argument, if given, must be ""--integral""\n");
            return 1;
        }
        integral = true;
        sscanf(argv[9], "%lld", &cn);
        sscanf(argv[10], "%lld", &cd);
    }

    // allocate the interpolation
    Poly* c = new Poly[nPieces];

    // determine the interpolation
    for (int i=0; i<nPieces; ++i)
    {
        BigFloat substart(i);
        substart.Mult(range).Add(start);
        BigFloat subfinish(substart);
        subfinish.Add(range);
        Chebyshev ch(F, NULL, substart, subfinish, nTerms*2);
        c[i].Copy(ch.GetPoly());
    }

    if (integral)
    {
        // integrate the interpolation to get the interpolation of the integral
        // C for the integral causes the low end of the range to evaluate to 0
        for (int i=0; i<nPieces; ++i)
        {
            for (int j=nTerms-1; j--;)
            {
                BigFloat x(c[i].At(j));
                x.Mult(range).Div(2).Div(j+1);
                c[i].SetAt(j+1, x);
            }
            c[i].SetAt(0, BigFloat(0));
        }
        BigFloat oldCoef(cn);
        oldCoef.Div(cd);
        for (int i=0; i<nPieces; ++i)
        {
            BigFloat coef(c[i].Eval(BigFloat(-1)));
            coef.Negate();
            BigFloat coef2(c[i].Eval(BigFloat(1)));
            coef.Add(oldCoef);
            if (i > 0)
            {
                coef.Add(c[i-1].At(0));
            }
            c[i].SetAt(0, coef);
            oldCoef.Copy(coef2);
        }
    }

    // print out functions implementing the interpolation
    printf("// %s %d/%d %d/%d %d %d %s %d/%d\n", name, sn, sd, fn, fd, nPieces, nTerms, integral ? "t" : "f", cn, cd);
    printf("static bool isInit%s = false;\n", name);
    printf("static const int nTerms%s = %d;\n", name, nTerms);
    printf("static double* coef%s = (double*)0;\n", name);
    printf("static double start%s;\n", name);
    printf("static double end%s;\n", name);
    printf("static double range%s;\n", name);
    printf("\n");
    printf("void Init%s()\n", name);
    printf("{\n");
    printf("    start%s = %d / static_cast<double>(%d);\n", name, sn, sd);
    printf("    range%s = %d / static_cast<double>(%d);\n", name, (fn*sd - sn*fd), (fd*sd*nPieces));
    printf("    coef%s = new double[%d];\n", name, nPieces*nTerms);
    for (int iPiece=0; iPiece<nPieces; ++iPiece)
    {
        int scale = 1; // we want to evaluate [-0.5, 0.5] but the coefficients are for [-1, 1]
        for (int iTerm=0; iTerm<nTerms; ++iTerm)
        {
            printf("    coef%s[%d] = ", name, iPiece*nTerms + iTerm);
            BigFloat ci(c[iPiece].At(iTerm));
            ci.Mult(scale).PrintDouble();
            printf(";\n");
            scale *= 2;
        }
        printf("\n");
    }
    printf("    isInit%s = true;\n", name);
    printf("}\n");
    printf("\n");
    printf("double %s(double x)\n", name);
    printf("{\n");
    printf("    if (!isInit%s)\n", name);
    printf("        Init%s();\n", name);
    printf("    x = (x-start%s)/range%s;\n", name, name);
    printf("    int piece = static_cast<int>(x);\n");
    printf("    double *c = &coef%s[piece * %d];\n", name, nTerms);
    printf("    x -= piece + 0.5;\n");
    printf("    double y = c[%d];\n", nTerms-1);
    for (int i=nTerms-1; i--; )
    {
        printf("    y = y * x + c[%d];\n", i);
    }
    printf("    return y;\n");
    printf("}\n");
    printf("\n");
    printf("#ifdef TESTFUNC\n");
    printf("#include <stdio.h>\n");
    printf("#include <math.h>\n");
    printf("#include \"bigFloat.h\"\n");
    printf("int main(int argc, char** argv)\n");
    printf("{\n");
    printf("    double iterations = 137;\n");
    printf("    for (int i=0; i<iterations; ++i)\n");
    printf("    {\n");
    printf("        double x = i/(iterations-1);\n");
    printf("        double y = -cos(x);\n");
    printf("        double z = %s(x);\n", name);
    printf("        BigFloat b(i);\n");
    printf("        double z2 = b.Div(iterations-1).Cos().Negate().ToDouble();\n");
    printf("        printf(\"%%g  (%%g %%g)  (%%g %%g) (%%g)\\n\", x, y, (z2-y), z, (z2-z), (y-z));\n");
    printf("    }\n");
    printf("}\n");
    printf("#endif /* TESTFUNC */\n");

    return 0;
}