//  ============================================================
//
//  ec.h
//
//  Copyright 2002, Dennis Meilicke and Rene Peralta
//
//  ============================================================
//
//  Description:
//
//  Experimental implementation of an elliptic curve class.
//  I would be really surprised if this worked...  Use at your
//  own risk.
//
//  Should this really be part of NTTL?
//
//  ============================================================
#ifndef __nttl_ec__
#define __nttl_ec__


#include <nttl/traits.h>
#include <nttl/gcd.h>


//
//  The ecPoint class should probably have a flag
//  indicating whether or not all of the coordinates
//  should be calculated.
//


template< class T > class ecPoint;


template< class T >
class ecCurve
{
    //
    //      y^2 = x^3 + A * x + B
    //
    friend class ecPoint<T>;

private:
    T _a;
    T _b;
    T _m;

public:
    ecCurve( ) : _a(0), _b(0) { }
    ecCurve( const T& m ) : _a(0), _b(0), _m(m) { }
    ecCurve( const T& a, const T& b, const T& m ) : _a(a), _b(b), _m(m) { }
    ecCurve( const ecCurve<T>& c ) { _a = c._a; _b = c._b; _m = c._m; }

    ecCurve& Random( size_t digits )
    {
        nttlTraits<T> t;

        _a = t.Random( digits );
        _b = t.Random( digits );

        return *this;
    }

    ecCurve& Random( T mod )
    {
        nttlTraits<T> t;

        _a = t.RandomMod( mod );
        _b = t.RandomMod( mod );

        return *this;
    }

    ecCurve& Random( ecPoint<T> p, size_t digits )
    {
        nttlTraits<T> t;

        _a = t.Random( digits );

        _b = ( t.Square( p._y ) % _m )
           - ( ( t.Square( p._x ) % _m ) * p._x % _m )
           - ( _a * p._x );
        _b %= _m;

        return *this;
    }

    bool operator==( const ecCurve<T> &c ) const
    {
        return ( (_a == c._a) && (_b == c._b) && (_m == c._m) );
    }

    bool operator!=( const ecCurve<T> &c ) const { return !operator==(c); }

    bool OnCurve( ecPoint<T> p )
    {
        nttlTraits<T> t;

        T lhs = t.Square( p._y ) % _m;

        T rhs = ( ( ( t.Square( p._x ) % _m ) * p._x % _m )
              + ( _a * p._x ) + _b ) % _m;

        return ( lhs == rhs );
    }
};


template< class T >
class ecPoint
{
    friend class ecCurve<T>;

private:
    T _x;
    T _y;
    T _z;
    ecCurve< T > *_c;

    T _Inverse( const T& x, const T& n ) const
    {
        T gcd, u, v;
        gcd = ExtendedGCD( x, n, &u, &v );
        //if ( gcd > 1 )
            //send_factor( gcd, n );
            //  throw exception?  What?
        return u;
    }

public:
    ecPoint( ) : _x(0), _y(0), _z(1) { }
    ecPoint( const ecCurve< T >& c ) : _x(0), _y(0), _z(1) { _c = (ecCurve<T> *)&c; }
    ecPoint( const ecPoint< T >& p )
                    { _x = p._x; _y = p._y; _z = p._z; _c = p._c; }
    ecPoint( const ecCurve<T>&c, const T& x, const T& y )
           : _x(x), _y(y), _z(1) { _c = (ecCurve<T> *)&c; }
    ecPoint( const ecCurve<T>&c, const T& x, const T& y, const T& z )
           : _x(x), _y(y), _z(z) { _c = (ecCurve<T> *)&c; }

    ecPoint<T>& Assign( const T& x, const T& y, const T& c )
                { _x = x; _y = y; _z = 1; _c = (ecCurve<T> *)&c; return *this; }
    ecPoint<T>& Assign( const T& x, const T& y, const T& z, const T& c )
                { _x = x; _y = y; _z = z; _c = (ecCurve<T> *)&c; return *this; }

    ecPoint<T>& operator=( const ecPoint<T>& p )
                { _x = p._x; _y = p._y; _z = p._z; _c = p._c; return *this; }

    ecPoint<T> operator+( const ecPoint<T>& p ) const
    {
        nttlTraits<T> t;

        if( p._c != _c ) { cerr << "different curves!" << endl; exit( 1 ); }

        ecPoint<T> C;

        if( ( _x == p._x ) && ( _y == (-p._y) ) )
        {
            //  Identity
            C.Assign( 0, 1, 0 );
            return C;
        }

        T lambda;
        if( ( _x == p._x ) && ( _y == p._y ) )
            lambda = ( 3 * t.Square( this->_x ) + _c->_a )
                   * _Inverse( 2 * _y, _c->_m );
        else
            lambda = ( _y - p._y) * _Inverse( _x - p._x, _c->_m );
        lambda = lambda % _c->_m;

        C._x = (lambda * lambda - _x - p._x) % _c->_m;
        C._y = (_c->_m - ( lambda * ( C._x - _x ) + _y ) ) % _c->_m;
        C._z = 1;

        return C;
    }

    ecPoint<T>& operator+=( const ecPoint<T>& p )
    {
        nttlTraits<T> t;

        if( p._c != _c ) { cerr << "different curves!" << endl; exit( 1 ); }

        if( ( _x == p._x ) && ( _y == (-p._y) ) )
        {
            //  Identity
            Assign( 0, 1, 0 );
            return *this;
        }

        T lambda;
        if( ( _x == p._x ) && ( _y == p._y ) )
            lambda = ( 3 * t.Square( this->_x ) + _c->_a )
                   * _Inverse( 2 * _y, _c->_m );
        else
            lambda = ( _y - p._y) * _Inverse( _x - p._x, _c->_m );
        lambda = lambda % _c->_m;

        T x(_x);
        _x = (lambda * lambda - x - p._x) % _c->_m;
        _y = (_c->_m - ( lambda * ( _x - x ) + _y ) ) % _c->_m;
        _z = 1;

        return *this;
    }

    ecPoint<T> operator*( long e )
    //
    //  This function takes on an unusual form so that we
    //  don't have to worry about the identity element.
    //
    //  This function should be parameterized differently,
    //  so that the exponent (the parameter e) can be of
    //  any type.  Does that mean removing this function 
    //  from the class (and making it stand-alone)?
    //
    {
        nttlTraits<long> t;

        ecPoint<T> acc( *this );
        while( t.IsEven( e ) )
        {
            acc += acc;
            e >>= 1;
        }

        ecPoint<T> a( acc );
        acc += acc;
        e >>= 1;

        while( e > 0 )
        {
            if( t.IsOdd( e ) )
                a += acc;
            acc += acc;
            e >>= 1;
        }

        return a;
    }

    ecPoint<T>& operator*=( long e )
    //
    //  See the notes for operator*
    //
    {
        nttlTraits<T> t;

        ecPoint<T> acc( *this );
        while( t.IsEven( e ) )
        {
            acc += acc;
            e >>= 1;
        }

        operator=( acc );
        acc += acc;
        e >>= 1;

        while( e > 0 )
        {
            if( t.IsOdd( e ) )
                operator+=( acc );
            acc += acc;
            e >>= 1;
        }

        return *this;
    }

    bool operator==( const ecPoint<T> &p ) const
    {
        return ( (_x == p._x) && (_y == p._y) && (_z == p._z) && (_c == p._c) );
    }

    bool operator!=( const ecPoint<T> &p ) const { return !operator==( p ); }
};


#endif