//  ============================================================
//
//  modSqrtLehmer.h
//
//  Copyright 2002, Dennis Meilicke and Rene Peralta
//
//  ============================================================
//
//  Description:
//
//  Modular square root using Lehmer's algorithm.  See [Knuth2],
//  pp ????
//
//  ============================================================
#ifndef __nttlHeader_modSqrtLehmer__
#define __nttlHeader_modSqrtLehmer__


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


#error modSqrtLehmer.h: This function is still under construction


template< typename T >
T
SqrtLehmer( const T& X, const T& p )
{
  T S;
  switch( ModSqrtCmn( a, p, &S ) )
  {
    case enumModSqrtQNR:    S = 0;
    case enumModSqrtFound:  return S;
  }

  T qnr = ModSqrtQNR( 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;
        }
    }
    
    S = ( v * j ) % p;
    return S;
}


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


#ifndef __nttlAlgorithm_modSqrt__
#define __nttlAlgorithm_modSqrt__


template< class T >
inline
T
Sqrt( const T& x, const T& p )
{
  return SqrtLehmer( x, p );
}


#endif  //  __nttlAlgorithm_modSqrt__


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


#endif  //  __nttlHeader_modSqrtLehmer__