// figure out how much more efficient chebyshev nodes are than equally
// spaced nodes where you interpolate between the middle nodes

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


Poly& Poly::Zero()
{
    for (int i=0; i<c_max; ++i)
        _a[i].Copy(9999);
    _len = 0;
    return *this;
}


Poly& Poly::One()
{
    _len = 1;
    _a[0].Copy(1);
    return *this;
}


Poly& Poly::SetAt(int index, const BigFloat& d)
{
    ASSERT(index < c_max, "too big");
    if (index >= _len)
        for (; _len<=index; ++_len)
            _a[_len].Zero();
    _a[index].Copy(d);
    while (_len > 0 && _a[_len-1].IsZero())
        --_len;
    return *this;
}


Poly& Poly::Set(const BigFloat& v1, const BigFloat& v2)
{
    _len = 2;
    _a[0].Copy(v1);
    _a[1].Copy(v2);
    return *this;
}


const BigFloat& Poly::At(int index) const
{
    return index < _len ? _a[index] : BigFloat::ConstZero();
}


BigFloat Poly::Eval(const BigFloat& x) const
{
    BigFloat result(0);
    for (int i=_len; i--;)
    {
        result.Mult(x);
        result.Add(_a[i]);
    }
    return result;
}


double Poly::EvalDouble(const BigFloat& x) const
{
    double xd = x.ToDouble();
    double result = 0.0;
    for (int i=_len; i--;)
    {
        result *= xd;
        result += _a[i].ToDouble();
    }
    return result;
}


Poly& Poly::Copy(const Poly& p)
{
    _len = p._len;
    for (int i=0; i<_len; ++i)
            _a[i].Copy(p._a[i]);
    return *this;
}


Poly& Poly::Add(const Poly& p)
{
    if (_len < p._len)
        for (int i=_len; i<p._len; i++)
            _a[i].Zero();
    for (int i=0; i<p._len; i++)
        _a[i].Add(p._a[i]);
    return *this;
}


Poly& Poly::Mult(const Poly& p)
{
    Poly q(*this);
    if (_len == 0 || p._len == 0)
    {
        _len = 0;
    }
    else
    {
        _len = _len + p._len - 1;
        ASSERT(_len <= c_max, "too big");
        for (int i=0; i<_len; ++i)
            _a[i].Zero();
        for (int i=0; i<p._len; ++i)
            for (int j=0; j<q._len; ++j)
            {
                BigFloat z(p._a[i]);
                z.Mult(q._a[j]);
                    _a[i+j].Add(z);
            }
    }
    return *this;
}


Poly& Poly::Mult(const BigFloat& d)
{
    for (int i=0; i<_len; ++i)
        _a[i].Mult(d);
    return *this;
}


void Poly::Show() const
{
    printf("poly len %d: ", _len);
    for (int i=0; i<_len; ++i)
        printf(" %g", _a[i].ToDouble());
}


void Lagrange::Init(
    BigFloat (*f)(const BigFloat&, void* context),
    void *context,
    const BigFloat& start,
    const BigFloat& finish,
    int n)
{
    _f = f;
    _context = context;
    _start.Copy(start);
    _finish.Copy(finish);
    _n = n;
    
        // support point indexes
    _s = new BigFloat[n+1];
    DefineSupport();
    
    // support point values
    _v = new BigFloat[n+1];
    for (int i=0; i<=n; ++i)
    {
        BigFloat x(_finish);
        x.Add(_start);
        BigFloat y(_finish);
        y.Sub(_start);
        y.Mult(_s[i]);
        x.Add(y);
        x.Div(2);
        _v[i].Copy((*_f)(x, _context));
    }
    
    // orthogonal polynomial
    _o = new Poly[n+1];
    Ortho();
    
    // interpolating polynomial
    _p.Zero();
    for (int iPoly=0; iPoly<=_n; ++iPoly)
    {
        for (int iDegree=0; iDegree<=_n; ++iDegree)
        {
            BigFloat coef(_p.At(iDegree));
            coef.Add(_o[iPoly].At(iDegree));
            _p.SetAt(iDegree, coef);
        }
    }
    
    // remove irrelevant terms
    while (_p.Length() > 0)
    {
        BigFloat relevance(_p.At(0));
        relevance.Add(_p.At(_p.Length()-1));
        if (relevance.Compare(_p.At(0)) != 0)
            break;
        _p.SetAt(_p.Length()-1, BigFloat::ConstZero());
    }
}


Lagrange::~Lagrange()
{
    delete[] _s;
    delete[] _v;
    delete[] _o;
}

    
// Given n points (a[i],v[i]), set *(p[i]) to the ith Lagrange polynomial
void Lagrange::Ortho()
{
    BigFloat minusOne(-1);
    for (int i=0; i<=_n; i++)
    {
        _o[i].One();
        for (int j=0; j<=_n; j++)
        {
            if (i==j)
                continue;
            
            Poly q(_s[j], minusOne);
            _o[i].Mult(q);
        }
        BigFloat scale(_v[i]);
        _o[i].Mult(scale.Div(_o[i].Eval(_s[i])));
    }
}

    
// Interpolate f(x)
BigFloat Lagrange::Eval(const BigFloat& x) const
{
    // convert x (between _start and _finish) to x2 (between -1 and 1)
    BigFloat range(_finish);
    range.Sub(_start);
    BigFloat x2(x);
    x2.Sub(_start).Mult(2).Div(range).Sub(1);
    
    // interpolate the value for x2
    return _p.Eval(x2);
}


// Interpolate f(x)
double Lagrange::EvalDouble(const BigFloat& x) const
{
    // convert x (between _start and _finish) to x2 (between -1 and 1)
    BigFloat range(_finish);
    range.Sub(_start);
    BigFloat x2(x);
    x2.Sub(_start).Mult(2).Div(range).Sub(1);
    
    // interpolate the value for x2
    return _p.EvalDouble(x2);
}


// absolute value of difference between interpolation and f at x
BigFloat Lagrange::Delta(const BigFloat& x) const
{
    BigFloat delta(Eval(x));
    delta.Sub((*_f)(x, _context));
    if (delta.Compare(BigFloat::ConstZero()) < 0)
        delta.Negate();
    return delta;
}


// absolute value of difference between interpolation and f at x
double Lagrange::DeltaDouble(const BigFloat& x) const
{
    double delta = EvalDouble(x);
    delta -= (*_f)(x, _context).ToDouble();
    if (delta < 0)
        delta = -delta;
    return delta;
}


// sample several points in the interval and return the worst delta seen
BigFloat Lagrange::WorstDelta() const
{
    static const int c_deltas = 13;
    BigFloat q[c_deltas];
    BigFloat inRange(_finish);
    inRange.Sub(_start).Div(2).Add(_start);
    q[0].Copy(Delta(inRange));
    int k = 1;
    for (int i=3; i<=7; i+=2)
    {
        for (int j=1; j<i; ++j)
        {
            inRange.Copy(_finish);
            inRange.Sub(_start).Mult(j).Div(i).Add(_start);
            q[k+j-1].Copy(Delta(inRange));
        }
        k += i-1;
    }
    ASSERT(k == c_deltas, "k=%d, c_deltas=%d\n", k, c_deltas);
    BigFloat worstDelta(q[0]);  // left-to-right, not pemdas $$$ got this far
    for (int i=1; i<c_deltas; ++i)
        if (q[i].Compare(worstDelta) > 0)
            worstDelta.Copy(q[i]);
    return worstDelta;
}



// sample several points in the interval and return the worst delta seen
double Lagrange::WorstDeltaDouble() const
{
    static const int c_deltas = 13;
    double q[c_deltas];
    BigFloat inRange(_finish);
    inRange.Sub(_start).Div(2).Add(_start);
    q[0] = DeltaDouble(inRange);
    int k = 1;
    for (int i=3; i<=7; i+=2)
    {
        for (int j=1; j<i; ++j)
        {
            inRange.Copy(_finish);
            inRange.Sub(_start).Mult(j).Div(i).Add(_start);
            q[k+j-1] = DeltaDouble(inRange);
        }
        k += i-1;
    }
    ASSERT(k == c_deltas, "k=%d, c_deltas=%d\n", k, c_deltas);
    double worstDelta = q[0];
    for (int i=1; i<c_deltas; ++i)
        if (q[i] > worstDelta)
            worstDelta = q[i];
    return worstDelta;
}




#ifdef TEST_BIGLAGRANGE

void Driver(const BigFloat& start, const BigFloat& finish, BigFloat (*f)(const BigFloat&, void*), const char* name)
{
    int n;
    void* context = NULL;
    BigFloat cWorst(0);
    BigFloat eWorst(0);
    for (n=1; n<30; ++n)
    {
        EqualSpacing e(f, context, start, finish, n);
        eWorst.Copy(e.WorstDelta());
        Chebyshev c(f, context, start, finish, n);
        cWorst.Copy(c.WorstDelta());
        printf("Degree %2d %s, Chebyshev=%g, EqualSpacing=%g ", n, name, cWorst.ToDouble(), eWorst.ToDouble());
        const Poly& p = c.GetPoly();
        BigFloat relevant(p.At(0));
        relevant.Add(p.At(n));
        if (relevant.Compare(p.At(0)) == 0)
            printf("c pointless ");
        const Poly& q = e.GetPoly();
        relevant.Copy(q.At(0)).Add(q.At(n));
        if (relevant.Compare(q.At(0)) == 0)
            printf("degree %d", q.Length()-1);
        else
        {
            BigFloat third(finish);
            third.Sub(start).Div(3).Add(start);
            EqualSpacing h(f, context, start, third, n);
            BigFloat hWorst(h.WorstDelta());
            printf(", ThirdSpacing=%g", hWorst.ToDouble());
        }
        printf("\n");
    }
}


void DriverDouble(const BigFloat& start, const BigFloat& finish, BigFloat (*f)(const BigFloat&, void*), const char* name)
{
    int n;
    void* context = NULL;
    double cWorst(0);
    double eWorst(0);
    for (n=1; n<15; ++n)
    {
        EqualSpacing e(f, context, start, finish, n);
        eWorst = e.WorstDeltaDouble();
        Chebyshev c(f, context, start, finish, n);
        cWorst = c.WorstDeltaDouble();
        printf("Degree %2d %s, Chebyshev=%g, EqualSpacing=%g ", n, name, log10(cWorst), log10(eWorst));
        const Poly& p = c.GetPoly();
        BigFloat relevant(p.At(0));
        relevant.Add(p.At(n));
        if (relevant.Compare(p.At(0)) == 0)
            printf("c pointless ");
        const Poly& q = e.GetPoly();
        relevant.Copy(q.At(0)).Add(q.At(n));
        if (relevant.Compare(q.At(0)) == 0)
            printf("degree %d", q.Length()-1);
        else
        {
            BigFloat third(finish);
            third.Sub(start).Div(3).Add(start);
            EqualSpacing h(f, context, start, third, n);
            double hWorst(h.WorstDeltaDouble());
            printf(", ThirdSpacing=%g", log10(hWorst));
        }
        printf("\n");
    }
}


BigFloat Ranga(const BigFloat& x, void*)
{
    BigFloat r(x);
    r.Mult(r);
    return x.IsNegative() ? r : r.Negate();
}


BigFloat Runge(const BigFloat& x, void*)
{
    BigFloat r(x);
    return r.Mult(r).Add(25).Invert();
}

BigFloat Sine(const BigFloat& x, void*)
{
    BigFloat r(x);
    return r.Sin();
}

BigFloat Exp(const BigFloat& x, void*)
{
    BigFloat r(x);
    return r.Exp();
}

int main()
{
    BigFloat a(-1);
    a.Div(100);
    BigFloat b(1);
    b.Div(100);
    try
    {
        DriverDouble(a, b, &Exp, "exp");
        DriverDouble(a, b, &Sine, "sine");
        DriverDouble(a, b, &Runge, "Runge");
    }
    catch( const std::exception & ex )
    {
        fprintf(stderr, "%s\n", ex.what());
    }
}

#endif // TEST_BIGLAGRANGE