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

#include "bigfloat.h"


// polynomials, with addition and multiplication and evaluation
// _a has _len terms, for x^^0 .. x^^(_len-1)
class Poly
{
public:
    Poly() { Zero(); }

    Poly(const Poly& p) { Copy(p); }

    Poly(const BigFloat& v1, const BigFloat& v2) { Set(v1, v2); }

    ~Poly() {}

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

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

    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 (_a[_len-1].IsZero() && _len > 0)
            --_len;
        return *this;
    }

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

    int Length() const
    {
        return _len;
    }

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

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

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

    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& 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& Mult(const BigFloat& d)
    {
        for (int i=0; i<_len; ++i)
            _a[i].Mult(d);
        return *this;
    }

private:
    static const int c_max = 40;
    BigFloat _a[c_max];
    int _len;
    
};



// An nth degree Lagrange interpolation for f from start to finish
class Lagrange
{
public:
    void Init(BigFloat (*f)(const BigFloat&), const BigFloat& start, const BigFloat& finish, int n)
    {
        _f = f;
        _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));
        }

        // 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());
        }
    }

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

    virtual void DefineSupport() = 0;
    
    // Given n points (a[i],v[i]), set *(p[i]) to the ith Lagrange polynomial
    void 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 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);
    }

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


    // sample several points in the interval and return the worst delta seen
    BigFloat 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].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;
    }


    // get the actual interpolation polynomial
    const Poly& GetPoly() const
    {
        return _p;
    }

    
protected:
    BigFloat (*_f)(const BigFloat&);  // the function to interpolate
    BigFloat  _start; // start of interval
    BigFloat  _finish; // end of interval
    int     _n;  // degree of interpolation polynomial
    BigFloat *_s;  // indexes of n+1 support points
    BigFloat *_v;  // values of n+1 support points
    Poly   *_o;  // orthogonal polynomials for n+1 support points
    Poly    _p;  // the Lagrange polynomial for this interval
};


// _n+1 Chebyshev nodes between -1 and 1
class Chebyshev : public Lagrange
{
public:
    Chebyshev(BigFloat (*f)(const BigFloat&), const BigFloat& start, const BigFloat& finish, int n)
    {
        Init(f, start, finish, n);
    }

    ~Chebyshev() {}

    void DefineSupport()
    {
        for (int i=0; i<=_n; ++i)
        {
            BigFloat pii(i);
            pii.Mult(BigFloat::Pi()).Div(_n).Cos().Negate();
            _s[i].Copy(pii);
        }
    }
};


// Equally spaced nodes, two apart, centered on 0
class EqualSpacing : public Lagrange
{
public:
    EqualSpacing(BigFloat (*f)(const BigFloat&), const BigFloat& start, const BigFloat& finish, int n)
    {
        Init(f, start, finish, n);
    }

    ~EqualSpacing() {}
    
// generate _n+1 equally-spaced nodes, each 2 apart, centered on 0
    void DefineSupport()
    {
        for (int i=0; i<=_n; ++i)
        {
            BigFloat si(i);
            si.Mult(2).Sub(_n);
            _s[i].Copy(si);
        }
    }
};



void Driver(const BigFloat& start, const BigFloat& finish, BigFloat (*f)(const BigFloat&), const char* name)
{
    int n;
    BigFloat cWorst(0);
    BigFloat eWorst(0);
    for (n=1; n<30; ++n)
    {
        EqualSpacing e(f, start, finish, n);
        eWorst.Copy(e.WorstDelta());
        Chebyshev c(f, 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, start, third, n);
            BigFloat hWorst(h.WorstDelta());
            printf(", ThirdSpacing=%g", hWorst.ToDouble());
        }
        printf("\n");
    }
}

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


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

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

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

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