#pragma once
#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();

    Poly& One();

    Poly& SetAt(int index, const BigFloat& d);

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

    int Length() const { return _len; }

    const BigFloat& At(int index) const;

    BigFloat Eval(const BigFloat& x) const;

    double EvalDouble(const BigFloat& x) const;

    Poly& Copy(const Poly& p);

    Poly& Add(const Poly& p);

    Poly& Mult(const Poly& p);

    Poly& Mult(const BigFloat& d);

    void Show() const;

private:
    static const int c_max = 60;
    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&, void* context),
        void* context,
        const BigFloat& start,
        const BigFloat& finish,
        int n);

    virtual ~Lagrange();

    virtual void DefineSupport() = 0;
    
    // Given n points (a[i],v[i]), set *(p[i]) to the ith Lagrange polynomial
    void Ortho();
    
    // Interpolate f(x)
    BigFloat Eval(const BigFloat& x) const;
    
    // Interpolate f(x)
    double EvalDouble(const BigFloat& x) const;
    
    // absolute value of difference between interpolation and f at x
    BigFloat Delta(const BigFloat& x) const;

    // absolute value of difference between interpolation and f at x
    double DeltaDouble(const BigFloat& x) const;

    // sample several points in the interval and return the worst delta seen
    BigFloat WorstDelta() const;

    // sample several points in the interval and return the worst delta seen
    double WorstDeltaDouble() const;

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

    const Poly& GetOrtho(int i) const { return _o[i]; }
    
protected:
    BigFloat (*_f)(const BigFloat&, void*);  // the function to interpolate
    void* _context;  // context for the function
    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&, void* context), void* context, const BigFloat& start, const BigFloat& finish, int n)
    {
        Init(f, context, 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&, void* context),
        void *context,
        const BigFloat& start,
        const BigFloat& finish,
        int n)
    {
        Init(f, context, 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);
        }
    }
};


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

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