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

#define EXTRA 5
#define SAMPLE 3

// G(t) = (1-cos(t))*(1-z*cos(t))^(-3/2), where z = 1 - e^^-y and y is context
BigFloat G(const BigFloat& t, void* context)
{
    
    BigFloat z(*(BigFloat*)context);
    z.Negate().Exp().Negate().Add(1);

    BigFloat b(t);
    b.Cos();
    BigFloat numerator(1);
    // The next line makes the numerator 1-cos, comment out the next line and the numerator is just 1.
    // numerator.Sub(b);
    b.Mult(z).Negate().Add(1);        
    BigFloat denominator(b);
    denominator.Mult(b).Mult(b).Sqrt();
    numerator.Div(denominator);

    return numerator;
}

// F(x) = f(y), where x = -ln(1-y) and f(y) = integral t=0..2pi of (1/2pi)*(1-cos(t))*(1-y*cos(t))^(-3/2) dt
BigFloat F(const BigFloat& x, void* context)
{
    BigFloat start(0);
    BigFloat finish(BigFloat::Pi());
    s8 nTerms = 30;

    const int minLogPieces = 5;
    const int maxLogPieces = 63;
    int logPieces = minLogPieces;
    BigFloat range(finish);
    range.Sub(start).Div(1<<logPieces);

    // accuracy desired for interpolation
    BigFloat small(2);
    small.Power(BigFloat(-64));

    // determine the interpolation
    BigFloat subStart(start);
    BigFloat startCoef(0);
    while (subStart.Compare(finish) < 0)
    {
        // interpolation polynomial for the current piece
        Poly c;
        bool lastPiece = false;
        BigFloat subFinish(subStart);
        subFinish.Add(range);
        for (;;)
        {
            Chebyshev ch(G, (void*)&x, subStart, subFinish, nTerms+EXTRA);
            c.Copy(ch.GetPoly());

            BigFloat biggestFirst(0);
            BigFloat biggestMiddle(0);
            BigFloat biggestLast(0);
            for (int i=0; i<SAMPLE; ++i)
            {
                if (biggestFirst.CompareAbsolute(c.At(i)) < 0)
                    biggestFirst.Copy(c.At(i));
                if (biggestMiddle.CompareAbsolute(c.At(i+nTerms/2)) < 0)
                    biggestMiddle.Copy(c.At(i+nTerms/2));
                if (biggestLast.CompareAbsolute(c.At(i+nTerms)) < 0)
                    biggestLast.Copy(c.At(i+nTerms));
            }
            biggestFirst.Mult(small);

            // too much error?
            BigFloat t(subStart);
            t.Add(range);

            // too little error?
            if (lastPiece)
                ;
            else if (biggestFirst.CompareAbsolute(biggestMiddle) > 0 && logPieces > minLogPieces)
            {
                BigFloat newSubFinish(subFinish);
                newSubFinish.Add(range);
                if (newSubFinish.Compare(finish) <= 0)
                {
                    logPieces--;
                    subFinish.Copy(newSubFinish);
                    range.Mult(2);
                    continue;
                }
            }
            else if (biggestFirst.CompareAbsolute(biggestLast) < 0 && logPieces < maxLogPieces)
            {
                logPieces++;
                range.Div(2);
                subFinish.Copy(subStart);
                subFinish.Add(range);
                continue;
            }
            else if (t.Compare(finish) > 0)
            {
                lastPiece = true;
                range.Copy(finish);
                range.Sub(subStart);
                subFinish.Copy(finish);
                continue;
            }

            break;
        }

        // integrate this piece
        for (int j=c.Length(); j-- > 0;)
        {
            BigFloat x(c.At(j));
            x.Mult(range).Div(2).Div(j+1);
            c.SetAt(j+1, x);
        }
        c.SetAt(0, BigFloat(0));

        // divide by 2pi
        c.Mult(BigFloat(2).Mult(BigFloat::Pi()).Invert());

        BigFloat coef(c.Eval(BigFloat(-1)));
        coef.Negate();
        BigFloat coef2(c.Eval(BigFloat(1)));
        startCoef.Add(coef);
        startCoef.Add(coef2);

        subStart.Copy(subFinish);
    }

    // evaluate the integral at 2pi (it will be twice its value at pi)
    BigFloat r(startCoef); // result
    r.Mult(2);
    
    return r;
}

int main(int argc, char** argv)
{
    const char* name = "Th1";
    static const int f = 40;  // finish
    static const int maxTerms = 20;
    BigFloat start(0);
    BigFloat finish(f);
    s8 nPieces[f];

    // number of pieces for each power of e
    nPieces[0] = 16;
    nPieces[1] = 8;
    nPieces[2] = 4;
    nPieces[3] = 4;
    for (int i=4; i<f; ++i)
        nPieces[i] = 2;

    int* nTerms[f];
    for (int i=0; i<f; ++i)
        nTerms[i] = new int[nPieces[i]];
    
    // accuracy desired for interpolation
    BigFloat small(2);
    small.Power(BigFloat(-56));

    // determine the interpolation
    printf("#include <math.h>\n");
    printf("static double **coef%s[%d];\n", name, f);
    printf("static double allCoef%s[] = {\n", name);
    for (int iInt=0; iInt<f; ++iInt)
    {
        BigFloat range(1);
        range.Div(nPieces[iInt]);
        BigFloat walk(iInt);
        BigFloat terms[maxTerms];
        for (int iPiece=0; iPiece<nPieces[iInt]; ++iPiece)
        {
            BigFloat subfinish(walk);
            subfinish.Add(range);
            Chebyshev ch(F, NULL, walk, subfinish, maxTerms+EXTRA);
            const Poly& p = ch.GetPoly();

            for (int iTerm=0; iTerm<maxTerms; ++iTerm)
            {
                BigFloat pi(p.At(iTerm));
                terms[iTerm].Copy(pi);
            }

            // figure out how big the early terms are
            BigFloat big(0);
            for (int i=0; i<maxTerms; ++i)
                if (big.CompareAbsolute(terms[i]) < 0)
                    big.Copy(terms[i]);
            big.Mult(small);

            // throw away the irrelevant tail
            nTerms[iInt][iPiece] = maxTerms;
            while (nTerms[iInt][iPiece]--)
            {
                if (big.CompareAbsolute(terms[nTerms[iInt][iPiece]]) < 0)
                {
                    nTerms[iInt][iPiece]++;
                    break;
                }
            }

            int scale = 1; // we want to evaluate [-0.5, 0.5] but the coefficients are for [-1, 1]
            for (int iTerm=0; iTerm<nTerms[iInt][iPiece]; ++iTerm)
            {
                terms[iTerm].Mult(scale);
                scale *= 2;
            }

            for (int iTerm = 0;  iTerm < nTerms[iInt][iPiece];  ++iTerm)
                printf("    %.18g,\n", terms[iTerm].ToDouble());
            printf("    // piece[%d][%d]  [%g,%g)  %d terms\n\n",
                   iInt, iPiece, walk.ToDouble(), subfinish.ToDouble(), nTerms[iInt][iPiece]);

            walk.Copy(subfinish);
        }
        printf("\n");
    }
    printf("};  // coef%s\n\n", name);

    printf("static int *nTerms%s[%d];\n", name, f);
    printf("static int allTerms%s[] = {\n", name);
    for (int iInt = 0;  iInt < f;  ++iInt)
    {
        for (int iPiece = 0;  iPiece < nPieces[iInt];  ++iPiece)
        {
            printf("    %d,\n", nTerms[iInt][iPiece]);
        }
        printf("\n");
    }
    printf("}; // allTerms%s\n\n", name);
    
    printf("static int nPieces%s[%d] = {\n", name, f);
    for (int i=0; i<f; ++i)
        printf("    %d,\n", nPieces[i]);
    printf("};  // nPieces%s\n\n", name);
    
    // print out function implementing the interpolation
    printf("static bool isInit%s = false;\n", name);
    printf("static int finish%s = %d;\n", name, f);
    printf("\n");
    printf("void Init%s()\n", name);
    printf("{\n");
    printf("    for (int iInt = 0, walkTerm = 0, walkCoef = 0;\n");
    printf("         iInt < finish%s;\n", name);
    printf("         walkTerm += nPieces%s[iInt], ++iInt)\n", name);
    printf("    {\n");
    printf("        nTerms%s[iInt] = &allTerms%s[walkTerm];\n", name, name);
    printf("        coef%s[iInt] = new double*[nPieces%s[iInt]];\n", name, name);
    printf("        for (int iPiece = 0;  iPiece < nPieces%s[iInt];  walkCoef += nTerms%s[iInt][iPiece], ++iPiece)\n",
           name, name);
    printf("        {\n");
    printf("            coef%s[iInt][iPiece] = &allCoef%s[walkCoef];\n", name, name);
    printf("        }\n");
    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\n", name);
    printf("    // Transform the domain to approach pole at 1 slower.\n");
    printf("    if (x < 0)\n");
    printf("        return nan(\"\");\n");
    printf("    x = -log(1-x);\n");
    printf("    if (x < 0 || x >= finish%s)\n", name);
    printf("        return nan(\"\");\n");
    printf("    int iInt = static_cast<int>(x);\n");
    printf("    x -= iInt;\n\n");
    printf("    // The 0.5 adjusts the range to [-0.5,0.5) instead of [0,1), so c[0] is the middle value.\n");
    printf("    // The number of pieces per int was adjusted manually.\n");
    printf("    x *= nPieces%s[iInt];\n", name);
    printf("    int iPiece = static_cast<int>(x);\n");
    printf("    x = x - iPiece - 0.5;\n\n");
    printf("    // start with the tiny terms and build up, to avoid rounding errors\n");
    printf("    double *c = coef%s[iInt][iPiece];\n", name);
    printf("    int iTerm = nTerms%s[iInt][iPiece]-1;\n", name);
    printf("    double y = c[iTerm];\n");
    printf("    while (iTerm--)\n");
    printf("        y = y * x + c[iTerm];\n");
    printf("    return y;\n");
    printf("}\n");
    printf("\n");
    printf("#ifdef TESTFUNC\n");
    printf("#include <stdio.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 z = %s(x);\n", name);
    printf("        printf(\"%%.18g => %%.18g\\n\", x, z);\n");
    printf("    }\n");
    printf("}\n");
    printf("#endif /* TESTFUNC */\n");

    return 0;
}