#include "biglagrange.h"

struct Point
{
    void Init(const BigFloat& x, const BigFloat& y, const BigFloat& z)
    {
        _x.Copy(x);
        _y.Copy(y);
        _z.Copy(z);
    }

    void Copy(Point p)
    {
        Init(p._x, p._y, p._z);
    }


    void Sub(Point p)
    {
        _x.Sub(p._x);
        _y.Sub(p._y);
        _z.Sub(p._z);
    }

    BigFloat Distance(Point p)
    {
        BigFloat d;
        BigFloat t;
        d.Add(t.Copy(_x).Sub(p._x).Mult(t));
        d.Add(t.Copy(_y).Sub(p._y).Mult(t));
        d.Add(t.Copy(_z).Sub(p._z).Mult(t));
        return d.Sqrt();
    }

    void Show()
    {
        printf("%g %g %g", _x.ToDouble(), _y.ToDouble(), _z.ToDouble());
    }
    
    BigFloat _x;
    BigFloat _y;
    BigFloat _z;
};


class Ring
{
public:
    // build a ring with radius 1, around the origin, in the (x,z) plane, with nMoons moons in it
    Ring(s8 nMoons, bool odd)
    {
        _nMoons = nMoons;
        _odd = odd;
        _ring = new Point[_nMoons];
        // make all the points for a ring of 
        for (s8 i=0; i<nMoons; i++)
        {
            BigFloat x(2*i+odd);
            x.Mult(BigFloatCache::_twoPi);
            x.Div(2*nMoons);
            BigFloat z(x);
            z.Sin();
            x.Cos();
            _ring[i].Init(x, 0, z);
        }
    }

    ~Ring()
    {
        delete[] _ring;
    }

    // find the force exerted by a unit ring of unit mass on a point exactly by nbody calculations
    Point NBodyForce(const Point& p) const
    {
        // find the acceleration it exerts on the point of reference
        Point force;
        for (s8 i=0; i<_nMoons; ++i)
        {
            // do not measure the attraction of a point to itself
            if (_ring[i]._x.Compare(p._x) == 0 &&
                _ring[i]._y.Compare(p._y) == 0 &&
                _ring[i]._z.Compare(p._z) == 0)
                continue;
            
            Point q(_ring[i]);
            q.Sub(p);
            BigFloat t;
            BigFloat d;
            d.Add(t.Copy(q._x).Mult(t));
            d.Add(t.Copy(q._y).Mult(t));
            d.Add(t.Copy(q._z).Mult(t));
            t.Copy(d);
            t.Mult(d.Sqrt());
            t.Mult(_nMoons);
            force._x.Add(q._x.Div(t));
            force._y.Add(q._y.Div(t));
            force._z.Add(q._z.Div(t));
        }
        return force;
    }

    
private:
    s8 _nMoons;
    Point* _ring;
    bool _odd;
};


struct RingContext
{
    const Ring* r;
    const Point* p;
};


static BigFloat Fx(const BigFloat& x)
{
    BigFloat fx(x);
    fx.Sub(1);
    fx.Invert();
    fx.Add(1);
    return fx;
}

static BigFloat Fy(const BigFloat& x, const BigFloat& y)
{
    BigFloat fy(y);
    BigFloat fx(x);
    fx.Sub(1);
    fx.Invert();
    fy.Mult(fx);
    return fy;
}

static BigFloat XNBodyX(const BigFloat& x, void* context)
{
    const Ring* r = ((RingContext*)context)->r;
    Point p(*((RingContext*)context)->p);
    p._y.Copy(Fy(x,p._y));
    p._x.Copy(Fx(x));
    Point accelOdd = r->NBodyForce(p);
    return accelOdd._x;
}

static BigFloat XNBodyY(const BigFloat& x, void* context)
{
    const Ring* r = ((RingContext*)context)->r;
    Point p(*((RingContext*)context)->p);
    p._y.Copy(Fy(x,p._y));
    p._x.Copy(Fx(x));
    Point accelOdd = r->NBodyForce(p);
    return accelOdd._y;
}

static BigFloat YNBodyX(const BigFloat& y, void* context)
{
    const Ring* r = ((RingContext*)context)->r;
    Point p(*((RingContext*)context)->p);
    p._y.Copy(Fy(p._x,y));
    p._x.Copy(Fx(p._x));
    Point accelOdd = r->NBodyForce(p);
    return accelOdd._x;
}

static BigFloat YNBodyY(const BigFloat& y, void* context)
{
    const Ring* r = ((RingContext*)context)->r;
    Point p(*((RingContext*)context)->p);
    p._y.Copy(Fy(p._x,y));
    p._x.Copy(Fx(p._x));
    Point accelOdd = r->NBodyForce(p);
    return accelOdd._y;
}

void Driver(s8 nMoons)
{
    // the point being acted upon
    Point p;
    BigFloat x(2);
    BigFloat y(-1);
    p.Init(x, y, BigFloat::ConstZero());
    Point p2(p);
    p2._x.Copy(Fx(x));
    p2._y.Copy(Fy(x, y));
    
    // make all the points for a ring of nMoons moons
    Ring oddRing(nMoons, true);
    Ring evenRing(nMoons, false);

    // find the acceleration it exerts on the point of reference
    Point accelEven = oddRing.NBodyForce(p2);
    Point accelOdd = evenRing.NBodyForce(p2);

    // report the result
    accelOdd.Show();
    printf("  diff %g\n", accelEven.Distance(accelOdd).ToDouble());

    RingContext ctx = {&oddRing, &p};
    BigFloat gap(1);
    gap.Div(256);
    BigFloat start(p._y);
    start.Sub(gap);
    BigFloat finish(p._y);
    finish.Add(gap);
    for (int i=1; i<20; ++i)
    {
        EqualSpacing interpolate(YNBodyY, &ctx, start, finish, i);
        interpolate.GetPoly().Show();
        printf("\n");
        BigFloat guess(interpolate.Eval(p._y));
        BigFloat exact(oddRing.NBodyForce(p2)._y);
        BigFloat diff(exact);
        diff.Sub(guess);
        p2.Show();
        printf(" degree %d, guess=%g, exact=%g, diff %g\n", i, guess.ToDouble(), exact.ToDouble(), diff.ToDouble());
    }
}

int main(int argc, char** argv)
{
    BigFloatCache::Init();
    
    ASSERT(argc == 2);
    s8 nMoons;
    sscanf(argv[1], "%ld", &nMoons);
    printf("nmoons: %d\n", nMoons);

    try
    {
        Driver(nMoons);
    }
    catch( const std::exception & ex )
    {
        fprintf(stderr, "%s\n", ex.what());
    }
    return 0;
}