//  ============================================================
//
//  jacobi.cpp
//
//  Copyright 2002, Dennis Meilicke and Rene Peralta
//
//  ============================================================
//
//  Description:
//
//  Implementation of the Jacobi function. [bressoud]
//
//  Note that this file is only temporary.  This code should be
//  moved to NTTL asap.
//
//  ============================================================
#include <lnv3/lnv3.h>


static ln
FactorPowerOfTwo( ln A, digit_t *e )
//
//  Description:
//
//  Compute K s.t. A = K * 2^e ;  Return K
//
//  This should be using the NTTL version...
//
{
    if( A == 0 )
    {
        (*e) = 0;
        return A;
    }

    for( (*e)=0 ; A.IsEven( ) ; (*e)++ )
        A >>= 1;

    return A;
}


// ---------------------------------------------------------------------- //


inline short
smallJacobi( const ln &n, const ln &modulus )
//
//  Description
//
//  Compute the Jacobi of n mod modulus, where n <= 2.
//  This is just a helper function for the main Jacobi
//  function, and is not intended to be exposed to the
//  outside world.
//
//  Warning
//
//  This function will fail for n > 2.
//  
{
  digit_t temp;
  digit_t digit = n.GetDigit( 0 );
  
  if( digit == 0 )
    return 0;

  if( digit == 1 )
    return 1;
  
  //
  //  If we are here, n == 2
  //
  //  The "& 7", below, is just a fast way to do "mod 8".
  //
  temp = modulus.GetDigit( 0 ) & 7;
  if( temp == 1 || temp == 7 )
    return 1;
  
  return -1;
}


// ---------------------------------------------------------------------- //


short
ln::Jacobi( const ln &Modulus ) const
//
//  Description
//
//  Calculate the value of the Jacobi symbol (*this / n).
//
//  Notes
//
//  Valid for all *this and odd n.
//
//  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]
//                
{
  ln Mod;
  ln k;
  ln a;
  digit_t s;    // 's' is returned by the 'GetEvenPower' function. 
  
  short J = 1;          // initial value for the Jacobi symbol

  Mod = Modulus;
  a = (*this) % Mod;
  
  while (1)
  {
    // boundary case
    if ( a <= 2 ) 
    {
      J = J * smallJacobi( a, Mod );
      return J;
    }

    if( a.IsOdd( ) )
    {
      //
      //  '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.
      //
      k = a;
      a = Mod % a;
      if( a == 1 ) 
         return J;
      Mod = k;
    }
    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.
      //  
      a = FactorPowerOfTwo( a, &s );
      if( s & 1 )  // if s is odd
      {
        //
        //  The "& 0x7", below, is just a fast way to do "mod 8".
        //
        digit_t tmp1 = Mod.GetDigit( 0 ) & 0x7;
        if( (tmp1 == 3) || (tmp1 == 5) )
          J = -J;
      }
    } 
  }
}