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

#include "base.h"
#include <math.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(double v1, double v2) { Set(v1, v2); }

    ~Poly() {}

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

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

    Poly& SetAt(int index, double d)
    {
        ASSERT(index < c_max, "too big");
        if (index >= _len)
            for (; _len<=index; ++_len)
                _a[_len] = 0;
        _a[index] = d;
        while (_a[_len-1] == 0.0 && _len > 0)
            --_len;
        return *this;
    }

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

    int Length() const
    {
        return _len;
    }

    double At(int index) const
    {
        return index < _len ? _a[index] : 0.0;
    }

    double Eval(double x) const
    {
        double result = 0.0;
        for (int i=_len; i--;)
        {
            result *= x;
            result += _a[i];
        }
        return result;
    }

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

    Poly& Add(const Poly& p)
    {
        if (_len < p._len)
            for (int i=_len; i<p._len; i++)
                _a[i] = 0;
        for (int i=0; i<p._len; i++)
            _a[i] += 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] = 0.0;
            for (int i=0; i<p._len; ++i)
                for (int j=0; j<q._len; ++j)
                    _a[i+j] += p._a[i] * q._a[j];
        }
        return *this;
    }

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

private:
    static const int c_max = 20;
    double _a[c_max];
    int _len;
    
};



// An nth degree Lagrange interpolation for f from start to finish
class Lagrange
{
public:
    void Init(double (*f)(double), double start, double finish, int n)
    {
        _f = f;
        _start = start;
        _finish = finish;
        _n = n;

        // support point indexes
        _s = new double[n+1];
        DefineSupport();

        // support point values
        _v = new double[n+1];
        for (int i=0; i<=n; ++i)
            _v[i] = (*_f)((_finish+_start)/2 + (_finish-_start)*_s[i]/2);

        // 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)
            {
                double coef = _p.At(iDegree);
                coef += _o[iPoly].At(iDegree);
                _p.SetAt(iDegree, coef);
            }
        }
        while (_p.Length() > 0 && _p.At(_p.Length()-1) + _p.At(0) == _p.At(0))
            _p.SetAt(_p.Length()-1, 0.0);
    }

    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()
    {
        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], -1.0);
                _o[i].Mult(q);
            }
            _o[i].Mult(_v[i]/_o[i].Eval(_s[i]));
        }
    }

    
    // Interpolate f(x)
    double Eval(double x) const
    {
        // convert x (between _start and _finish) to x2 (between -1 and 1)
        double x2 = 2.0 * (x - _start)/(_finish - _start) - 1.0;

        // interpolate the value for x2
        return _p.Eval(x2);
    }

    
    // absolute value of difference between interpolation and f at x
    double Delta(double x) const
    {
        double delta = Eval(x) - (*_f)(x);
        if (delta < 0.0)
            delta = -delta;
        return delta;
    }


    // sample several points in the interval and return the worst delta seen
    double WorstDelta() const
    {
        static const int c_deltas = 13;
        double q[c_deltas];
        q[0] = Delta(_start+(_finish-_start)/2);
        int k = 1;
        for (int i=3; i<=7; i+=2)
        {
            for (int j=1; j<i; ++j)
                q[k+j] = Delta(_start+j*(_finish-_start)/i);
            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;
    }


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

    
protected:
    double (*_f)(double);  // the function to interpolate
    double  _start; // start of interval
    double  _finish; // end of interval
    int     _n;  // degree of interpolation polynomial
    double *_s;  // indexes of n+1 support points
    double *_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(double (*f)(double), double start, double finish, int n)
    {
        Init(f, start, finish, n);
    }

    ~Chebyshev() {}

    void DefineSupport()
    {
        for (int i=0; i<=_n; ++i)
            _s[i] = -cos(i*3.1415926535897932384626433/_n);
        
    }
};


// Equally spaced nodes, two apart, centered on 0
class EqualSpacing : public Lagrange
{
public:
    EqualSpacing(double (*f)(double), double start, double 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)
            _s[i] = -_n + 2*i;
    }
};


double EquivalentSpacing(
    double start,
    double finish,
    int n,
    double (*f)(double))
{
    Chebyshev* c = new Chebyshev(f, start, finish, n);
    EqualSpacing* e = new EqualSpacing(f, start, finish, n);
    double cWorst = c->WorstDelta();
    double eWorst = e->WorstDelta();
    if (eWorst == 0.0)
        return 0.0;
    else if (cWorst == 0.0)
        return 1000.0;
    else if (eWorst > cWorst)
    {
        // get within half
        double prev = 2.0;
        double next = 1.0;

        // binary search
        double delta = prev - next;
        for (int i=0; i<20; ++i)
        {
            prev = next;
            delta /= 2;
            if (eWorst > cWorst)
                next -= delta;
            else
                next += delta;
            delete e;
            e = new EqualSpacing(f, start, start + (finish-start)*next, n);
            eWorst = e->WorstDelta();
        }

        return next;
    }

}


void Driver(double start, double finish, double (*f)(double), const char* name)
{
    int n;
    double cWorst = 0.0;
    double eWorst = 0.0;
    for (n=1; n<15; ++n)
    {
        EqualSpacing e(f, start, finish, n);
        eWorst = e.WorstDelta();
        Chebyshev c(f, start, finish, n);
        cWorst = c.WorstDelta();
        printf("Degree %2d %s, Chebyshev=%g, EqualSpacing=%g ", n, name, log10(cWorst), log10(eWorst));
        const Poly& p = c.GetPoly();
        if (p.At(n) + p.At(0) == p.At(0))
            printf("c pointless ");
        const Poly& q = e.GetPoly();
        if (q.At(n) + q.At(0) == q.At(0))
            printf("degree %d", q.Length()-1);
        else
        {
            EqualSpacing h(f, start, start+(finish-start)/3, n);
            double hWorst = h.WorstDelta();
            printf(", ThirdSpacing=%g", log10(hWorst));
        }
        printf("\n");
    }
}

double Runge(double x)
{
    return 1.0/(25.0 + x*x);
}

int main()
{
    try
    {
        Driver(-0.01, 0.01, &exp, "exp  ");
        Driver(-0.01, 0.01, &sin, "sine ");
        Driver(-0.01, 0.01, &Runge, "Runge");
    }
    catch( const std::exception & ex )
    {
        fprintf(stderr, "%s\n", ex.what());
    }
}