//  ============================================================
//
//  modSqrt.h
//
//  Copyright 2002, Dennis Meilicke and Rene Peralta
//
//  ============================================================
//
//  Description:
//
//  NTTL Implementation of:
//      LehmerSqrt
//      Sqrt
//
//  Note:
//
//  ModSqrtRecursive should not be called from user code.
//
//  ============================================================
#ifndef __nttl_modSqrt__
#define __nttl_modSqrt__


#include <nttl/traits.h>
#include <nttl/inverse.h>
#include <nttl/factorPowerOfTwo.h>
#include <nttl/jacobi.h>

#error modSqrtOld.h: This header does not work.

#ifdef UNDER_CONSTRUCTION

template< class T, class TMOD >
TMOD
LehmerSqrt( const T& X, const TMOD& p )
//
//  Description:  Calculate modular square root of X mod p.
//  Algorithm:    Lehmer's algorithm (citation?)
//
{
    if( Jacobi( X, p ) == -1 )
    {
        cerr << "<ERROR-0012>[nt1]. Quadratic nonresidue in SqrtMod."
            << endl;
        exit(0);
    }
    
    TMOD x = X % p;
    if (x == 0) return 0;
    if (x == 1) return 1;  // this takes care of p == 2
    
    if ( (p % 4) == 3 ) return Power( x, (p+1)/4, p );
    if ( (p % 8) == 5 )
    {
        TMOD aux = Power( x, (p-1)/4, p ); // this is a square root of 1 mod p
        
        if ( aux == 1 ) return Power( x, (p+3)/8 , p );

        return ( ((p + 1) / 2) * Power( 4*x, (p+3)/8 , p ) ) % p;
    }
    
    /* now p ==  1 mod 8 */
    
    static LargeNumber j,y,w,m;
    long k;
    static LargeNumber x4;
    
    
    /* find h, the smallest number such that h*h - 4*x is 
    a quadratic non-residue mod p */
    
    x4 = 4*x; // for efficiency
    for ( h = 1 ; (h.Square() - x4).Jacobi(p) == 1 ; ++h);
    
    j = (p + 1)/2;
    
    /* initialize */
    
    m = x;
    v = h;
    w = (h.Square() - 2*m) % p;
    
    /* main loop */
    
    for (k = j.NumBits() - 2 ; k >= 0 ; k--)
    {
        aux = 2*m;
        
        y = (v * w - h * m) % p;
        v = (v.Square() - aux ) % p;
        w = (w.Square() -  x * aux ) % p;
        m = ( m.Square() ) % p;
        
        if ( j.Bit(k) == 0 )  w = y;
        else
        {
            v = y;
            m = ( x * m ) % p;
        }
    }
    
    return ( v * j ) % p;
}

#endif


//  ============================================================


//
//  Is there any way to make this "static"?  That is, un-callable
//  from the outside world?
//


template< class T >
T
PowerOf2( T* p, long e )
{
  T ret = 1;
  if( e )
    ret = ( *p << (e - 1) ) ;
  return ret;
}


template< class T >
T
ModSqrtRecursive( T x, const T& p, const T& qnr)
//
//  Description
//
//  Recursive helper function for the modular square root 
//  function
//  Notes
//
//  This should be converted to a non-recursive version.  Since
//  this is tail-recursive, this should be trivial.
//
//  The parameterization of the template should be changed so
//  that 'p' (the modulus) can be of a different type.
//
{
    //cout << "recursize call" << endl;
    //cout << "    x = " << x << endl;
    //cout << "    p = " << p << endl;
    //cout << "    q = " << qnr << endl;

    nttlTraits<T> t;

    long a;
    //cout << "before FactorPowerOfTwo" << endl;
    T s = FactorPowerOfTwo( ( p - 1 ), &a );
    //cout << "after FactorPowerOfTwo" << endl;
    
    //cout << "before h" << endl;
    //cout << "    x = " << x << endl;
    //cout << "    s = " << s << endl;
    //cout << "    p = " << p << endl;
    T h = t.Power( x, s, p );
    //cout << "after h.  h = " << h << endl;
    if( h == 1 )
        return t.Power( x, (s+1)/2, p );
    
    long k = 0;
    do
    {
        h = t.Square( h ) % p;
        //cout << "h = " << h << endl;
        k++;
    } while( h != 1 );
    //cout << "after k" << endl;
    
    T v;
    PowerOf2( &v, a - k - 1 );
    
    T n1 = t.Power( qnr, v, p );
    T n2 = t.Square( n1 ) % p;
    n1   = Inverse( n1, p );
    
    return ( n1 * ModSqrtRecursive( (x*n2) % p, p, qnr ) % p );
}


//  ============================================================


//
//  Can we do this so that the type of the modulus is different
//  than the type of the number?  That would be very nice...
//


template< class T >
T
Sqrt( const T &X, const T& p )
//
//  Description
//
//  Calculate a modular square root.
//
//  Algorithm
//  ============================================================
//
//  My notes say that this algorithm is due to Adleman, Manders, 
//  and Miller.  I can't find a citation for this.  It bears
//  some resemblance to the Tonelli and Shanks algorithm in 
//  [COHEN] (pg33), but I'm not sure.  Either way, we should
//  look at the algorithm in [COHEN], and consider implementing
//  it.
//
//  To do
//
//  This method uses a recursive helper function.  This should
//  at least be changed to be non-recursive.  It would be nice
//  to get rid of the helper function completely.
//                
{
    nttlTraits< T > t;

    if( X.Jacobi( p ) == -1 )
    {
        //  QNR. Throw an exception?
        cerr << "nttl:Sqrt: Quadratic nonresidue in SqrtMod." << endl;
        cerr << "X = " << X << endl;
        cerr << "p = " << p << endl;
        exit(0);
    }
    
    T x = X % p;
    if( x == 0 ) return 0;
    if( x == 1 ) return 1; // Note that this takes are of p == 2
    
    if( ( p % 4 ) == 3 ) return t.Power( x, (p+1)/4, p );
    if( ( p % 8 ) == 5 )
    {
        T aux = t.Power( x, (p-1)/4, p ); // this is a square root of 1 mod p
        
        if( aux == 1 ) return t.Power( x, (p+3)/8, p );

        return ( ((p + 1) / 2) * t.Power( (4*x), (p+3)/8, p ) ) % p;
    }
    
    //
    //  If we are here, then p == 1 mod 8
    //
    //  Find a Quadratic Nonresidue (QNR) mod p
    //
    T qnr;
    for( qnr = 2 ; qnr.Jacobi( p ) == 1 ; ++qnr ) ;
    
    return ModSqrtRecursive( x, p, qnr );
}


#endif