//  ============================================================
//
//  modSqrtTonelli.h
//
//  Copyright 2002, Dennis Meilicke and Rene Peralta
//
//  ============================================================
//
//  Description:
//
//  See [Menezes,van Oorschot,Vanstone], pg 100, algorithm 3.34
//
//  [Menezes,van Oorschot,Vanstone] says this algorithm is
//  derived from Koblitz, and mentioned in Bach and Shallit.
//
//  ============================================================
#ifndef __nttlHeader_modSqrtTonelli__
#define __nttlHeader_modSqrtTonelli__


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


template< typename T >
T
SqrtTonelli( T a, const T& p )
{
  nttlTraits<T> tr;

  if( a > p ) a %= p;
  T pMinus1 = p - 1;

  T r;
  switch( ModSqrtCmn( a, p, &r ) )
  {
    case enumModSqrtQNR:       r = 0;
    case enumModSqrtFound:     return r;
    case enumModSqrtNotFound:  break;
  }

  T b = ModSqrtQNR( p );

  long s;
  T t = FactorPowerOfTwo( pMinus1, &s );

  T aInv = Inverse( a, p );

  T c = tr.Power( b, t, p );

  r = tr.Power( a, (t+1)/2, p );

  T dExp(1);
  if( s > 0 ) dExp <<= (s - 1);

  for( long i=1; i<s ; i++ )
  {
    dExp >>= 1;
    T d = tr.Power( r * r * aInv, dExp, p );
    if( d == pMinus1 )  //  == -1 mod p
      r = r * c % p;
    c = c * c % p;
  }

  return r;
}


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


#ifndef __nttlAlgorithm_modSqrt__
#define __nttlAlgorithm_modSqrt__


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


#endif  //  __nttlAlgorithm_modSqrt__


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


#endif  //  __nttlHeader_modSqrtTonelli__