//  ============================================================
//
//  gcdEuclid.h
//
//  Copyright 2002, Dennis Meilicke and Rene Peralta
//
//  ============================================================
//
//  Description:
//
//  Implementation of the Euclidean GCD algorithm.
//
//  ============================================================
//
//  Functions:
//      GCDEuclid
//      EGCDEuclid
//
//  ============================================================
#ifndef __nttlHeader_gcdEuclid__
#define __nttlHeader_gcdEuclid__


#include <nttl/traits.h>


template< class T >
T
GCDEuclid( const T& a, const T& b ) 
//
//  Description
//
//  Return the Greatest Common Denominator of two numbers.
//
//  Notes
//
//  This could be sped up quite a bit by swapping address,
//  rather than copying LargeNumbers.
//
//  This should be replaced by Lehmer's routine.
// 
{
  T g0, g1, g2;
  
  g0 = ( b > 0 ) ? b : -b;
  g1 = ( a > 0 ) ? a : -a;
  
  if ( g1 > g0 )
  {
    T swap = g0;
    g0 = g1;
    g1 = swap;
  }
  
  while( g1 != 0 )
  {
    g2 = g0 % g1;
    g0 = g1;
    g1 = g2;
  }
  
  return g0;
}


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


#ifndef __nttlAlgorithm_gcd__
#define __nttlAlgorithm_gcd__


template< class T >
T
GCD( const T& a, const T& b )
{
  return GCDEuclid( a, b );
}


#endif  //  __nttlAlgorithm_gcd__


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


template< class T >
T
EGCDEuclid( const T& a, const T& b, T *u, T *v )
//
//  Description
//
//  Determine u, v, d such that au + bv = d
//
//  Algorithm
//
//  From [COHEN], p 16.
//
//  Notes
//
//  This implementation is straight from the book, but it seems
//  that it does too many "expensive" calculations (div, mod,
//  mult).
//
//  This function has been "patched" to work with negative 
//  numbers.  These patches should be reviewed, to see if there
//  is a more efficient solution.
//
{
  nttlTraits<T> t;

  *u = 1;
  T d = t.Abs( a );

  if( b == 0 )
  {
	*v = 0;
	return d;
  }

  T v1 = 0;
  T v3 = t.Abs( b );

  T t1, t3, q;

  while( v3 != 0 )
  {
	q = d / v3;
	t3 = d % v3;

	t1 = *u - q * v1;
	*u = v1;
	d  = v3;
	v1 = t1;
	v3 = t3;
  }

  *v = ( d - t.Abs(a) * *u ) / t.Abs(b);
  if (a < 0) *u = -(*u);
  if (b < 0) *v = -(*v);
  return d;
}


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


#ifndef __nttlAlgorithm_egcd__
#define __nttlAlgorithm_egcd__


template< class T >
T
EGCD( const T& u, const T& v, T* a, T* b )
{
  return EGCDEuclid( u, v, a, b );
}


#endif  //  __nttlAlgorithm_egcd__


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


#endif  //   __nttlHeader_gcdEuclid__