//  ============================================================
//
//  gcdRSB.h
//
//  Copyright 2002, Dennis Meilicke and Rene Peralta
//
//  ============================================================
//
//  Description:
//
//  Implementation of the right-shift binary GCD algorithm.
//
//  ============================================================
//
//  Functions:
//      GCDRSB - Right-shift binary GCD algorithm
//
//  ============================================================
//
//  References:
//
//  Jonathan Sorenson.  Two Fast GCD Algorithms.  Journal of
//  Algorithms 16, 110-144 (1994).
//
//  ============================================================
//
//  Acknowledgements:
//
//  Thanks to Fethiye Akblut for work on this algorithm.
//
//  ============================================================
//
//  To Do:
//
//  This algorithm could really benefit from an efficient
//  implementation of the factorPowerOfTwo algorithm.  To that
//  end, factorPowerOfTwo should be made a traits method.  The
//  non-traits version of the function (in factorPowerOfTwo.h)
//  would simply call the traits method.
//
//  For ln3, a first order -= function would speed this up, as
//  well.  Right now, -= is second order, implemented in terms
//  of regular -.
//
//  In the return expression, we only shift if g!=0.  What does
//  the standard say about x<<0 for a POD x?  How is it actually
//  implemented?  ln3 should do the same thing, and this 
//  algorithm should assume the standard behaviour.
//
//  ============================================================


#ifndef __nttlHeader_gcdRSB__
#define __nttlHeader_gcdRSB__


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


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


template< class T >
T
GCDRSB( T u, T v )
{
  nttlTraits<T> t;

  long gu;
  for( gu=0 ; t.IsEven( u ) ; gu++ )
    u >>= 1;

  long gv;
  for( gv=0 ; t.IsEven( v ) ; gv++ )
    v >>= 1;

  long g = ( gu > gv ) ? gv : gu;

  while( 1 )
  {
    if( u >= v )
    {
      u -= v;
      if( u == 0 ) return (g==0) ? v : (v<<=g);
      while( t.IsEven( u ) ) u >>= 1;
    }
    else
    {
      v -= u;
      while( t.IsEven( v ) ) v >>= 1;
    }
  }
}


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


#ifndef __nttlAlgorithm_gcd__
#define __nttlAlgorithm_gcd__


template< class T >
T
GCD( const T& u, const T& v )
{
  return GCDRSB( u, v );
}


#endif  //   __nttlAlgorithm_gcd__


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




template< class T >
T
EGCDRSBbad( T u, T v, T* u1, T* u2 )
{
  nttlTraits<T> t;

  long gu;
  for( gu=0 ; t.IsEven( u ) ; gu++ )
    u >>= 1;

  long gv;
  for( gv=0 ; t.IsEven( v ) ; gv++ )
    v >>= 1;

  long g = ( gu > gv ) ? gv : gu;

  *u1=1; *u2=0;  T u3=u;
  T v1=v,  v2=(1-u), v3=v;
  T t1=0,  t2=-1,    t3=-v;

  while( 1 )
  {
    while( t.IsEven( t3 ) )
	{
	  if( t.IsOdd( t1 ) || t.IsOdd( t2 ) )
	  {
		  t1 += v;  t2 -= u;
	  }
	  t2 >>= 1; t2 >>= 1; t3 >>= 1;
	}

	if( t3 > 0 )
	{  *u1 = t1; *u2 = t2; u3 = t3; }
	else
	{  v1 = v - t1;  v2 = (-u) - t2;  v3 = -t3;  }

	t1 = *u1 - v1;  t2 = *u2 - v2;  t3 = u3 - v3;

	if( t3 == 0 )
		return (g)?(u3 << g):u3;

	if( t1 < 0 )
	{
		t1 += v;  t2 -= u;
	}
  }
}


template< class T >
T
EGCDRSB( T u, T v, T* u1, T* u2 )
{
  nttlTraits<T> t;

  long gu;
  for( gu=0 ; t.IsEven( u ) ; gu++ )
    u >>= 1;

  long gv;
  for( gv=0 ; t.IsEven( v ) ; gv++ )
    v >>= 1;

  long g = ( gu > gv ) ? gv : gu;

  *u1=1; T u3=u;
  T v1=v,  v3=v;
  T t1=0,  t3=-v;

  while( 1 )
  {
    while( t.IsEven( t3 ) )
	{
	  if( t.IsOdd( t1 ) )
		  t1 += v;
	  t1 >>= 1; t3 >>= 1;
	}

	if( t3 > 0 )
	{  *u1 = t1; u3 = t3; }
	else
	{  v1 = v - t1;  v3 = -t3;  }

	t1 = *u1 - v1;  t3 = u3 - v3;

	if( t3 == 0 )
		return (g)?(u3 << g):u3;

	if( t1 < 0 )
		t1 += v;
  }
}


#endif  //  __nttlHeader_gcdRSB__