//  ============================================================
//
//  jacobi.h
//
//  Copyright 2002, Dennis Meilicke and Rene Peralta
//
//  ============================================================
//
//  Description:
//
//  NTTL Implementation of:
//      Jacobi
//
//  Note:
//
//  Jacobi2 should not be called from user code.
//
//  ============================================================
#ifndef __nttlHeader_jacobi__
#define __nttlHeader_jacobi__


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


template< class T >
short
Jacobi2( const T &n, const T &modulus )
//
//  Description: Single precision version of smallJacobi().
//  Warning:     This function is only valid for n <= 2.
//  
{
  if( n == 0 )
    return 0;

  if( n == 1 )
    return 1;

  //
  //  If we are here, n == 2
  //
  //T temp = modulus & 7;  // i.e., temp = modulus % (digitSigned_t) 8
  nttlTraits<T> t;
  short temp = t.Mod8( modulus );
  if( (temp == 1) || (temp == 7) )
    return 1;
  
  return -1;
}



template< class TNUM, class TMOD >
short
Jacobi( const TNUM &N, const TMOD &Modulus )
//
//  Description:  Calculate the value of the Jacobi symbol (*this / n).
//  Notes:        Valid for all N and odd Modulus.
//  Algorithm:    This is basically the "ordinary" algorithm.  We
//                should try implementing one of the faster
//                algorithms.
//
//                Makes use of the following identities (from [Flath]:
//
//                (-1/m) = (-1)^(m-1)/2 = /  1 if m==1 mod 4
//                                        \ -1 if m==3 mod 4
//
//                (2/m)  = (-1)^(m^2-1)/8 = /  1 if m==+-1 mod 8
//                                          \ -1 if m==+-3 mod 8
//
//                (n/m)  = (-1)^(m-1)(n-1)/2
//
//                       = /  (m/n) if m==1 or n==1 mod 4
//                         \ -(m/n) if m==n==3 mod 4
//                
//  References:   [K2], exercise 4.5.4.23.
//                [Shall]
//                [Flath]
//                
{
  nttlTraits<TMOD> t;

  //
  //  Initial value of the Jacobi symbol.  For return-value-optimization,
  //  this should be the only variable that is returned from this function.
  //
  short J = 1;

  //
  //  Temporary variable.
  //  This is declared here, so that its constructor is called
  //  only once.
  //
  TMOD prevA;
  
  TMOD Mod = Modulus;
  TMOD a = N % Mod;
    
  while (1)
  {
    //
    //  Base case.
    //
    if ( a <= 2 ) 
    {
      J = J * Jacobi2( a, Mod );
      return J;
    }

    if( t.IsOdd( a ) )
    {
      //
      //  'a' is odd
      //
      //  then (a/m) = (-1)^((a-1)(m-1)/4)(m mod a /a)
      //  so calculate (-1)^x
      //  Then swap a and m, and reduce m mod a.
      //
      //  Note the ((a-1)(m-1)/4) is odd only when a%4=3 and
      //  m%4=3.
      //  
      if( ( a % 4 == 3 ) && ( Mod % 4 == 3 ) )
        J = -J;
      
      //
      //  exchange a and n.  After reducing a, if a=1, save some
      //  time, and exit right away.
      //
      prevA = a;
      a = Mod % a;
      if( a == 1 ) 
        return J;
      Mod = prevA;
    }
    else
    {
      //
      //  'a' is even.
      //
      //  Since (pq/m) = (p/m)(q/m), factor powers of 2 out of
      //  'a'.  Then
      //
      //  Since (2/m)   = (-1)^((m^2 - 1)/8), then 
      //  If 's' is even, then J is unchanged.  If 's' is odd,
      //  then (2/m)^s = (2/m).
      //  (m^2 - 1)/8 is odd only if n%8 = 3 or 5.
      //  
      long s;
      a = FactorPowerOfTwo( a, &s );
      if( s & 1 )  // if s is odd
      {
        short temp = t.Mod8( Mod );
        if( (temp == 3) || (temp == 5) )
          J = -J;
      }
    } 
  }

  //
  //  Note that we will never get here.
  //  All returns are from within the while loop.
  //
}


#endif  //  __nttlHeader_jacobi__