#ifndef EUCLIDH
#define EUCLIDH

/* ************************************************************************ *
 *                                                                          *
 *                 Euclid.h  (C) 1999    Michael R. Toth                    *
 *                                                                          *
 *    Some templated functions for operations in Euclidean domains. The     *
 *  class parameter "Euclidean" should implement  +, -, *, /, %, +=, *=,    *
 *  /=, %=, =, ==, and != . It should be constructable from an "int".       *
 *                                                                          *
 * ************************************************************************ */         

template<class Euclidean> inline
bool Associates( const Euclidean& x, const Euclidean& y )
               { return  x==0? y==0: y==0? false: x%y==0? y%x==0: false; }

template<class Euclidean> inline
bool Congruent( const Euclidean& x, const Euclidean& y, const Euclidean& Chi )
     // True if  x=y (mod Chi) .
              { return  Chi==0? x==y: (x-y)%Chi==0; }

template<class Euclidean>
Euclidean GCD( Euclidean m, Euclidean n )
             { for ( Euclidean d; m!=0; n=d )  d=m, m=n%d;  return n; }

template<class Euclidean> inline
Euclidean LCM( const Euclidean& p, Euclidean q )
             { q/=GCD(p,q), q*=p;  return q; }

template<class Euclidean> inline
bool Unit( const Euclidean& x )  // True if  x  is invertable.
         { return  x==0? false: Euclidean(1)%x==0; }

template<class Euclidean>
Euclidean Inverse( Euclidean& u, Euclidean Chi =0 )
          // (Return value)*(old u) = G.C.D.{old u, Chi} = new u  (mod Chi) .
                 {
                   if ( Chi==0 )  return Euclidean(1);
                   u%=Chi;
                   if ( u==0 )  { u=Chi;  return Euclidean(0); }
                   Euclidean  S0=0, S1=1, s;
                   for ( Euclidean Rm=Chi%u; Rm!=0; Chi=u, u=Rm, Rm=Chi%u )
                         s=S1, S1=S0-Chi/u*s, S0=s;
                   return S1;
                 }

template<class Euclidean>
bool Invert( Euclidean& x, const Euclidean& Chi =0 )
     // Set  x:=x[inv] (mod Chi) .
           {
             if ( Chi==0 )
                  if ( Unit(x) )  { x=Euclidean(1)/x;  return true; }
                       else  return false;
             Euclidean  Gcd=x;
             x=Inverse(Gcd,Chi);
             if ( !Unit(Gcd) )  return false;
             x/=Gcd, x+=Chi, x%=Chi;
             return true;
           }

template<class Euclidean>
Euclidean GcdReplace( Euclidean& x, Euclidean& y )
          // (old x)*(new x)+(old y)*(new y) = G.C.D.{old x, old y} = return value.
                    {
                      Euclidean gcd=x;
                      Euclidean z=Inverse(gcd,y);
                      if ( y!=0 )  y=(gcd-x*z)/y;
                      x=z;
                      return gcd;
                    }

#endif EUCLIDH