synaps/arithm/BC.h

/*********************************************************************
*      This file is part of the source code of SYNAPS kernel.        *
*   Author(s): B. Mourrain, GALAAD, INRIA                            *
**********************************************************************
History: 
$Id: BC.h,v 1.1 2005/07/11 10:39:47 mourrain Exp $
**********************************************************************/
#ifndef SYNAPS_ARITHM__BC_H_INCLUDED_
#define SYNAPS_ARITHM__BC_H_INCLUDED_

#ifndef NDEBUG
        #include <cassert>
#endif

//========================== end of GMP custom allocation =====================
//#include <synaps/Basic.h>
//#include <synaps/gmp.h>

#ifdef win32
        #define const
#endif

#include <synaps/init.h>

__BEGIN_NAMESPACE_SYNAPS


template <class T> struct Scl<T>
#ifdef PL_CMM
: public GcObject
#endif
{
public:
  Scl<T> ( );
  Scl<T> (const Scl<T>& b);
  Scl<T> (signed long sl);
  Scl<T> (unsigned long ul);
  Scl<T> (int si);
  Scl<T> (const char* string, int base = 10);
  Scl<T> (  MP_INT new_z ) { z = new_z; }

  ~Scl<T> ();
  friend ostream&       operator << (ostream& os, const Scl<T>& b);
  friend istream&       operator >> (istream& is, Scl<T>& b);

  friend unsigned long  BigIntToUL (const Scl<T>& b);
  friend signed   long  BigIntToSL (const Scl<T>& b);

  friend size_t         log (const Scl<T>& b);
  friend size_t         log (const Scl<T>& b, int base);

  friend int            sign    (const Scl<T>& b);

  friend int            compare (const Scl<T>& b1, const Scl<T>& b2);
  friend int            compare (const Scl<T>& b, unsigned long ul);
  friend int            compare (const Scl<T>& b, long sl);

  friend bool           IsPositive      (const Scl<T>& b);
  friend bool           IsNegative      (const Scl<T>& b);
  friend bool           IsZero          (const Scl<T>& b);
  friend bool           IsOne           (const Scl<T>& b);
  friend bool           IsMinusOne      (const Scl<T>& b);
  friend bool           IsOdd           (const Scl<T>& b);
  friend bool           IsEven          (const Scl<T>& b);
  friend bool           IsPerfectSquare (const Scl<T>& b);
  friend bool           IsProbablyPrime (const Scl<T>& b, int reps = 25);

  bool                  operator == (const Scl<T>& rhs) const;
  bool                  operator != (const Scl<T>& rhs) const;
  bool                  operator >  (const Scl<T>& rhs) const;
  bool                  operator >= (const Scl<T>& rhs) const;
  bool                  operator <  (const Scl<T>& rhs) const;
  bool                  operator <= (const Scl<T>& rhs) const;

  bool                  operator == (long sl) const;
  bool                  operator != (long sl) const;
  bool                  operator >  (long sl) const;
  bool                  operator >= (long sl) const;
  bool                  operator <  (long sl) const;
  bool                  operator <= (long sl) const;

  bool                  operator == (int si) const;
  bool                  operator != (int si) const;
  bool                  operator >  (int si) const;
  bool                  operator >= (int si) const;
  bool                  operator <  (int si) const;
  bool                  operator <= (int si) const;

  bool                  operator == (unsigned long ul) const;
  bool                  operator != (unsigned long ul) const;
  bool                  operator >  (unsigned long ul) const;
  bool                  operator >= (unsigned long ul) const;
  bool                  operator <  (unsigned long ul) const;
  bool                  operator <= (unsigned long ul) const;

  //
  // -------- functions that modify at least one argument or `*this' ----------
  //

  void                  operator =  (const Scl<T>& rhs);
  void                  operator += (const Scl<T>& rhs);
  void                  operator -= (const Scl<T>& rhs);
  void                  operator *= (const Scl<T>& rhs);
  void                  operator /= (const Scl<T>& rhs);
  void                  operator %= (const Scl<T>& rhs);
  
  void                  operator =  (unsigned long rhs);
  void                  operator += (unsigned long rhs);
  void                  operator -= (unsigned long rhs);
  void                  operator *= (unsigned long rhs);
  void                  operator /= (unsigned long rhs);
  void                  operator %= (unsigned long rhs);

  void                  operator =  (long rhs);
  void                  operator += (long rhs);
  void                  operator -= (long rhs);
  void                  operator *= (long rhs);
  void                  operator /= (long rhs);
  void                  operator %= (long rhs);
  
  void                  operator ++ ( );
  void                  operator -- ( );

  void                  Div2Exp   (unsigned long exponent_of_2);
  void                  GCD       (const Scl<T>& b2);
  void                  Mod2Exp   (unsigned long exponent_of_2);
  void                  Mul2Exp   (unsigned long exponent_of_2);
  void                  PowMod    (const Scl<T>& exp, const Scl<T>& m);
  void                  PowMod    (unsigned long exp, const Scl<T>& m);
  void                  abs       ( );
  void                  factorial (unsigned long n);
  void                  negate    ( );
  void                  pow       (unsigned long exp);
  void                  pow       (unsigned long base, unsigned long exp);
  void                  quo       (const Scl<T>& divisor);
  void                  quo       (unsigned long divisor);
  void                  rem       (const Scl<T>& divisor);
  unsigned long         rem       (unsigned long divisor);
  void                  sqrt      ( );

  friend void           SqrtRem (Scl<T>& sqrt, Scl<T>& rem, const Scl<T>& b);
  
  friend void           QuoRem  (Scl<T>& q, Scl<T>& r,
                                 const Scl<T>& divdend, const Scl<T>& divisor);

  friend unsigned long  QuoRem  (Scl<T>& q, Scl<T>& r,
                                 const Scl<T>& divdend, unsigned long divisor);

  friend void           DivMod  (Scl<T>& q, Scl<T>& r,
                                 const Scl<T>& divdend, const Scl<T>& divisor);

  friend void           DivMod  (Scl<T>& q, Scl<T>& r, const Scl<T>& divdend,
                                 unsigned long divisor);

  friend void           ExtGCD  (Scl<T>& gcd, Scl<T>& a, Scl<T>& b,
                                 const Scl<T>& x, const Scl<T>& y);

  friend void           HalfExtGCD (Scl<T>& gcd, Scl<T>& a,
                                    const Scl<T>& x, const Scl<T>& y);


#ifdef PL_CMM
  void traverse();
#endif
  mpz z;
};


//
// -------- the corresponding `const' functions -----------------------------
//

Scl<T>&          operator - (const Scl<T>& b);

Scl<T>         operator + (const Scl<T>& b1, const Scl<T>& b2);
Scl<T>&        operator - (const Scl<T>& b1, const Scl<T>& b2);
Scl<T>&        operator * (const Scl<T>& b1, const Scl<T>& b2);
Scl<T>&        operator / (const Scl<T>& b1, const Scl<T>& b2);
Scl<T>&        operator % (const Scl<T>& b1, const Scl<T>& b2);

Scl<T>&        operator + (const Scl<T>& b, unsigned long ul);
Scl<T>&        operator - (const Scl<T>& b, unsigned long ul);
Scl<T>&        operator * (const Scl<T>& b, unsigned long ul);
Scl<T>&        operator / (const Scl<T>& b, unsigned long ul);
Scl<T>&        operator % (const Scl<T>& b, unsigned long ul);

Scl<T>&        operator + (unsigned long ul, const Scl<T>& b);
Scl<T>&        operator * (unsigned long ul, const Scl<T>& b);

Scl<T>&        operator + (const Scl<T>& b, long sl);
Scl<T>&        operator - (const Scl<T>& b, long sl);
Scl<T>&        operator * (const Scl<T>& b, long sl);
Scl<T>&        operator / (const Scl<T>& b, long sl);
Scl<T>&        operator % (const Scl<T>& b, long sl);

Scl<T>&        operator + (long sl, const Scl<T>& b);
Scl<T>&        operator * (long sl, const Scl<T>& b);

Scl<T>&        Div2Exp (const Scl<T>& b, unsigned long exponent_of_2);
Scl<T>&        GCD     (const Scl<T>& b1, const Scl<T>& b2);
Scl<T>&        gcd     (const Scl<T>& b1, const Scl<T>& b2);
Scl<T>&        Mod2Exp (const Scl<T>& b, unsigned long exponent_of_2);
Scl<T>&        Mul2Exp (const Scl<T>& b, unsigned long exponent_of_2);
Scl<T>&        PowMod  (const Scl<T>& b, const Scl<T>& exp, const Scl<T>& m);
Scl<T>&        PowMod  (const Scl<T>& b, unsigned long exp, const Scl<T>& m);
Scl<T>&        abs     (const Scl<T>& b);
// Scl<T>&        negate  (const Scl<T>& b);
Scl<T>&        pow     (const Scl<T>& base, unsigned long exp);
Scl<T>&        quo     (const Scl<T>& dividend, const Scl<T>& divisor);
Scl<T>&        quo     (const Scl<T>& dividend, unsigned long divisor);
Scl<T>&        rem     (const Scl<T>& dividend, const Scl<T>& divisor);
Scl<T>&        sqrt    (const Scl<T>& b);
unsigned long  rem     (const Scl<T>& b, unsigned long divisor);

Scl<T>&        BigIntFactorial (unsigned long n);


//-------------------------- Scl<T> (...) -------------------------------------

inline
Scl<T>::Scl<T> ( )
{
  mpz_init(&z);
}

inline
Scl<T>::Scl<T> (const Scl<T>& b)
{
  mpz_init_set(&z, &b.z);
}

inline
Scl<T>::Scl<T> (signed long int sl)
{
  mpz_init_set_si(&z, sl);
}

inline
Scl<T>::Scl<T> (unsigned long int ul)
{
  mpz_init_set_ui(&z, ul);
}

inline
Scl<T>::Scl<T> (int si)
{
  mpz_init_set_si(&z, (long) si);
}

inline
Scl<T>::Scl<T> (const char *string, int base)
{
  if (mpz_init_set_str(&z, string, base))
    cerr << "Scl<T>::Scl<T>: The string " << string 
         << " is not a valid number in base " << base << endl;
}


//-------------------------- Scl<T>::~Scl<T> () -------------------------------

inline
Scl<T>::~Scl<T> ( )
{
#if defined(PL_CMM) // || defined(PL_BARTLETT)
// PL_CMM || Bartlett
// It wouldn't be wrong to call `mpz_clear' but since it only calls
// `_PossoGMPDealloc', which does nothing for Bartlett and CMM, this
// is a small optimization.
// PL_CMM || Bartlett end
#else
  mpz_clear(&z);
#endif
}


//-------------------------- BigIntToUL ---------------------------------------

inline unsigned long
BigIntToUL (const Scl<T>& b)
{
  return mpz_get_ui(&b.z);
}


//-------------------------- BigIntToSL ---------------------------------------

inline signed long
BigIntToSL (const Scl<T>& b)
{
  return mpz_get_si(&b.z);
}


//-------------------------- log(...) -----------------------------------------

inline size_t
log (const Scl<T>& b)
{
  return mpz_size(&b.z);
}

inline size_t
log (const Scl<T>& b, int base)
{
  assert(base >= 2 && base <= 36);
  return mpz_sizeinbase(&b.z, base);
}

//-------------------------- sign ---------------------------------------------

inline int
sign (const Scl<T>& b)
{
  if (b.z._mp_size == 0)
    return 0;
  else
    return b.z._mp_size > 0 ? 1 : -1;
}


//-------------------------- compare (...) ------------------------------------

inline int
compare (const Scl<T>& b1, const Scl<T>& b2)
{
  return mpz_cmp(&b1.z, &b2.z);
}

inline int
compare (const Scl<T>& b, unsigned long ul)
{
  return mpz_cmp_ui(&b.z, ul);
}

inline int
compare (const Scl<T>& b, long sl)
{
  return mpz_cmp_si(&b.z, sl);
}


//-------------------------- IsPositive etc. ----------------------------------

inline bool
IsPositive (const Scl<T>& b)
{
  return b.z._mp_size > 0;
}

inline bool
IsNegative (const Scl<T>& b)
{
  return b.z._mp_size < 0;
}

inline bool
IsZero (const Scl<T>& b)
{
  return b.z._mp_size == 0;
}

inline bool
IsOne(const Scl<T>& b)
{
  return b.z._mp_size == 1 && b.z._mp_d[0] == 1;
}

inline bool
IsMinusOne (const Scl<T>& b)
{
  return b.z._mp_size == -1 && b.z._mp_d[0] == 1;
}

inline bool
IsOdd (const Scl<T>& b)
{
  return b.z._mp_size != 0 && (b.z._mp_d[0] & ((mp_limb_t) 1));
}

inline bool
IsEven (const Scl<T>& b)
{
  return b.z._mp_size == 0 || (b.z._mp_d[0] ^ ((mp_limb_t) 0));
}

inline bool
IsPerfectSquare (const Scl<T>& b)
{
  return mpz_perfect_square_p(&b.z);
}

inline bool
IsProbablyPrime (const Scl<T>& b, int reps /* = 25 */)
{
  return mpz_probab_prime_p(&b.z, reps);
}


//-------------------------- Scl<T> ==, !=, >, >=, <, <= ----------------------

inline bool
Scl<T>::operator == (const Scl<T>& rhs) const
{
  return mpz_cmp(&z, &rhs.z) == 0;
}

inline bool
Scl<T>::operator != (const Scl<T>& rhs) const
{
  return mpz_cmp(&z, &rhs.z) != 0;
}

inline bool
Scl<T>::operator > (const Scl<T>& rhs) const
{
  return mpz_cmp(&z, &rhs.z) > 0;
}

inline bool
Scl<T>::operator >= (const Scl<T>& rhs) const
{
  return mpz_cmp(&z, &rhs.z) >= 0;
}

inline bool
Scl<T>::operator < (const Scl<T>& rhs) const
{
  return mpz_cmp(&z, &rhs.z) < 0;
}

inline bool
Scl<T>::operator <= (const Scl<T>& rhs) const
{
  return mpz_cmp(&z, &rhs.z) <= 0;
}


//-------------------------- long ==, !=, >, >=, <, <= ------------------------

inline bool
Scl<T>::operator == (long sl) const
{
  return mpz_cmp_si(&z, sl) == 0;
}

inline bool
Scl<T>::operator != (long sl) const
{
  return mpz_cmp_si(&z, sl) != 0;
}

inline bool
Scl<T>::operator > (long sl) const
{
  return mpz_cmp_si(&z, sl) > 0;
}

inline bool
Scl<T>::operator >= (long sl) const
{
  return mpz_cmp_si(&z, sl) >= 0;
}

inline bool
Scl<T>::operator < (long sl) const
{
  return mpz_cmp_si(&z, sl) < 0;
}

inline bool
Scl<T>::operator <= (long sl) const
{
  return mpz_cmp_si(&z, sl) <= 0;
}

//-------------------------- int ==, !=, >, >=, <, <= ---------------


inline bool
Scl<T>::operator == (int si) const
{
  return mpz_cmp_si(&z, (long) si) == 0;
}

inline bool
Scl<T>::operator != (int si) const
{
  return mpz_cmp_si(&z, (long) si) != 0;
}

inline bool
Scl<T>::operator > (int si) const
{
  return mpz_cmp_si(&z, (long) si) > 0;
}

inline bool
Scl<T>::operator >= (int si) const
{
  return mpz_cmp_si(&z, (long) si) >= 0;
}

inline bool
Scl<T>::operator < (int si) const
{
  return mpz_cmp_si(&z, (long) si) < 0;
}

inline bool
Scl<T>::operator <= (int si) const
{
  return mpz_cmp_si(&z, (long) si) <= 0;
}


//-------------------------- unsigned long ==, !=, >, >=, <, <= ---------------


inline bool
Scl<T>::operator == (unsigned long ul) const
{
  return mpz_cmp_ui(&z, ul) == 0;
}

inline bool
Scl<T>::operator != (unsigned long ul) const
{
  return mpz_cmp_ui(&z, ul) != 0;
}

inline bool
Scl<T>::operator > (unsigned long ul) const
{
  return mpz_cmp_ui(&z, ul) > 0;
}

inline bool
Scl<T>::operator >= (unsigned long ul) const
{
  return mpz_cmp_ui(&z, ul) >= 0;
}

inline bool
Scl<T>::operator < (unsigned long ul) const
{
  return mpz_cmp_ui(&z, ul) < 0;
}

inline bool
Scl<T>::operator <= (unsigned long ul) const
{
  return mpz_cmp_ui(&z, ul) <= 0;
}

//-------------------------- Scl<T> =, +=, -=, *=, /=, %= ---------------------

inline void
Scl<T>::operator =  (const Scl<T>& rhs)
{
  if (this != &rhs)
    mpz_set(&z, &rhs.z);
}

inline void
Scl<T>::operator += (const Scl<T>& rhs)
{
  mpz_add(&z, &z, &rhs.z);
}

inline void
Scl<T>::operator -= (const Scl<T>& rhs)
{
  mpz_sub(&z, &z, &rhs.z);
}

inline void
Scl<T>::operator *= (const Scl<T>& rhs)
{
  mpz_mul(&z, &z, &rhs.z);
}

inline void
Scl<T>::operator /= (const Scl<T>& rhs)
{
//  mpz_div(&z, &z, &rhs.z);
  mpz_divexact(&z, &z, &rhs.z);
}

inline void
Scl<T>::operator %= (const Scl<T>& rhs)
{
  mpz_mod(&z, &z, &rhs.z);
}


//-------------------------- unsigned long =, +=, -=, *=, /=, %= --------------

inline void
Scl<T>::operator =  (unsigned long ul)
{
  mpz_set_ui(&z, ul);
}

inline void
Scl<T>::operator += (unsigned long ul)
{
  mpz_add_ui(&z, &z, ul);
}

inline void
Scl<T>::operator -= (unsigned long ul)
{
  mpz_sub_ui(&z, &z, ul);
}

inline void
Scl<T>::operator *= (unsigned long ul)
{
  mpz_mul_ui(&z, &z, ul);
}

inline void
Scl<T>::operator /= (unsigned long ul)
{
  mpz_div_ui(&z, &z, ul);
}

inline void
Scl<T>::operator %= (unsigned long ul)
{
  mpz_mod_ui(&z, &z, ul);
}


//-------------------------- long =, +=, -=, *=, /=, %= -----------------------

inline void
Scl<T>::operator =  (long sl)
{
  mpz_set_si(&z, sl);
}

inline void
Scl<T>::operator += (long sl)
{
  if (sl >= 0)
    mpz_add_ui(&z, &z, (unsigned long) sl);
  else
    mpz_sub_ui(&z, &z, ((unsigned long) -sl));
}

inline void
Scl<T>::operator -= (long sl)
{
  if (sl >= 0)
    mpz_sub_ui(&z, &z, (unsigned long) sl);
  else
    mpz_add_ui(&z, &z, (unsigned long) sl);
}

inline void
Scl<T>::operator *= (long sl)
{
  if (sl >= 0)
    mpz_mul_ui(&z, &z, (unsigned long) sl);
  else
    {
      z._mp_size = -z._mp_size;
      mpz_mul_ui(&z, &z, ((unsigned long) -sl));
    }
}

inline void
Scl<T>::operator /= (long sl)
{
  if (sl >= 0)
    mpz_div_ui(&z, &z, (unsigned long) sl);
  else
    {
      z._mp_size = -z._mp_size;
      mpz_div_ui(&z, &z, ((unsigned long) -sl));
    }
}

// inline void
// Scl<T>::operator %= (long sl)
// {
//   mpz_mod_si(&z, &z, sl);
// }


//-------------------------- prefix ++, and -- --------------------------------

inline void
Scl<T>::operator ++ ( )
{
  mpz_add_ui(&z, &z, 1);
}

inline void
Scl<T>::operator -- ( )
{
  mpz_sub_ui(&z, &z, 1);
}


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

inline void
Scl<T>::Div2Exp (unsigned long exponent_of_2)
{
  mpz_div_2exp(&z, &z, exponent_of_2);
}

inline void
Scl<T>::GCD (const Scl<T>& b2)
{
    mpz_gcd (&z, &z, &b2.z);
}


inline void
Scl<T>::Mod2Exp (unsigned long exponent_of_2)
{
  mpz_mod_2exp(&z, &z, exponent_of_2);
}

inline void
Scl<T>::Mul2Exp (unsigned long exponent_of_2)
{
  mpz_mul_2exp(&z, &z, exponent_of_2);
}

inline void
Scl<T>::PowMod (const Scl<T>& exp, const Scl<T>& m)
{
  mpz_powm(&z, &z, &exp.z, &m.z);
}

inline void
Scl<T>::PowMod (unsigned long exp, const Scl<T>& m)
{
  mpz_powm_ui(&z, &z, exp, &m.z);
}

inline void
Scl<T>::abs ( )
{
  if (z._mp_size < 0)
    z._mp_size = - z._mp_size;
}

inline void
Scl<T>::factorial (unsigned long n)
{
  mpz_fac_ui(&z, n);
}

inline void
Scl<T>::negate ( )
{
  z._mp_size = - z._mp_size;
}

inline void
Scl<T>::pow (unsigned long exp)
{
  mpz_pow_ui(&z, &z, exp);
}

inline void
Scl<T>::pow (unsigned long base, unsigned long exp)
{
  mpz_ui_pow_ui(&z, base, exp);
}

inline void
Scl<T>::quo (const Scl<T>& divisor)
{
  mpz_mdiv (&z, &z, &divisor.z);
}

inline void
Scl<T>::quo (unsigned long divisor)
{
  mpz_mdiv_ui(&z, &z, divisor);
}

inline void
Scl<T>::rem (const Scl<T>& divisor)
{
  mpz_mmod(&z, &z, &divisor.z);
}

inline unsigned long
Scl<T>::rem (unsigned long divisor)
{
  return mpz_mmod_ui(&z, &z, divisor);
}

inline void
Scl<T>::sqrt ( )
{
  mpz_sqrt(&z, &z);
}

inline void
SqrtRem (Scl<T>& sqrt, Scl<T>& rem, const Scl<T>& b)
{
  mpz_sqrtrem(&sqrt.z, &rem.z, &b.z);
}
 
inline void
QuoRem (Scl<T>& q, Scl<T>& r, const Scl<T>& divdend, const Scl<T>& divisor)
{
  mpz_mdivmod(&q.z, &r.z, &divdend.z, &divisor.z);
}

inline unsigned long
QuoRem (Scl<T>& q, Scl<T>& r, const Scl<T>& divdend, unsigned long divisor)
{
  return mpz_mdivmod_ui(&q.z, &r.z, &divdend.z, divisor);
}

inline void
DivMod (Scl<T>& q, Scl<T>& r, const Scl<T>& divdend, const Scl<T>& divisor)
{
  mpz_divmod(&q.z, &r.z, &divdend.z, &divisor.z);
}

inline void
DivMod (Scl<T>& q, Scl<T>& r, const Scl<T>& divdend, unsigned long divisor)
{
  mpz_divmod_ui(&q.z, &r.z, &divdend.z, divisor);
}

inline void
ExtGCD (Scl<T>& gcd, Scl<T>& a, Scl<T>& b, const Scl<T>& x, const Scl<T>& y)
{
  mpz_gcdext(&gcd.z, &a.z, &b.z, &x.z, &y.z);
}

inline void
HalfExtGCD (Scl<T>& gcd, Scl<T>& a, const Scl<T>& x, const Scl<T>& y)
{
    mpz_gcdext(&gcd.z, &a.z, 0, &x.z, &y.z);
}


//
// -------- the corresponding `const' functions -----------------------------
//

//#define NO_NRV 1
#if defined(__GNUG__) && !defined(NO_NRV)
#define NRVAL(name) return name
#define NRCODE(code) code
#define NONRCODE(code) 
#else
#define NRVAL(name) 
#define NRCODE(code)
#define NONRCODE(code) code
#endif


inline Scl<T>
operator - (const Scl<T>& b) NRVAL(r(b))
{
  NRCODE(r.negate();)
  NONRCODE(Scl<T> r(b); r.negate(); return r;)
}

inline Scl<T>
operator + (const Scl<T>& b1, const Scl<T>& b2) NRVAL(r(b1))
{
  NRCODE(r += b2;)
  NONRCODE(Scl<T> r(b1); r += b2; return r;)
}

inline Scl<T>
operator - (const Scl<T>& b1, const Scl<T>& b2) NRVAL(r(b1))
{
  NRCODE(r -= b2;)
  NONRCODE(Scl<T> r(b1); r -= b2; return r;)
}

inline Scl<T>
operator * (const Scl<T>& b1, const Scl<T>& b2) NRVAL(r(b1))
{
  NRCODE(r *= b2;)
  NONRCODE(Scl<T> r(b1); r *= b2; return r;)    
}

inline Scl<T>
operator / (const Scl<T>& b1, const Scl<T>& b2) NRVAL(r(b1))
{
  NRCODE(r /= b2;)
  NONRCODE(Scl<T> r(b1); r /= b2; return r;)
}

inline Scl<T>
operator % (const Scl<T>& b1, const Scl<T>& b2) NRVAL(r(b1))
{
  NRCODE(r %= b2;)
  NONRCODE(Scl<T> r(b1); r %= b2; return r;)
}

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

inline Scl<T>
operator + (const Scl<T>& b, unsigned long ul) NRVAL(r(b))
{
  NRCODE(r += ul;)
  NONRCODE(Scl<T> r(b); r += ul; return r;)
}

inline Scl<T>
operator - (const Scl<T>& b, unsigned long ul) NRVAL(r(b))
{
  NRCODE(r -= ul;)
  NONRCODE(Scl<T> r(b); r -= ul; return r;)
}

inline Scl<T>
operator * (const Scl<T>& b, unsigned long ul) NRVAL(r(b))
{
  NRCODE(r *= ul;)
  NONRCODE(Scl<T> r(b); r *= ul; return r;)
}

inline Scl<T>
operator / (const Scl<T>& b, unsigned long ul) NRVAL(r(b))
{
  NRCODE(r /= ul;)
  NONRCODE(Scl<T> r(b); r /= ul; return r;)
}

inline Scl<T>
operator % (const Scl<T>& b, unsigned long ul) NRVAL(r(b))
{
  NRCODE(r %= ul;)
  NONRCODE(Scl<T> r(b); r %= ul; return r;)
}

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

inline Scl<T>
operator + (unsigned long ul, const Scl<T>& b)
{
  return b + ul;
}

inline Scl<T>
operator * (unsigned long ul, const Scl<T>& b)
{
  return b * ul;
}

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

inline Scl<T>
operator + (const Scl<T>& b, long sl)  NRVAL(r(b))
{
  NRCODE(r += sl;)
  NONRCODE(Scl<T> r(b); r += sl; return r;)
}


inline Scl<T>
operator - (const Scl<T>& b, long sl)  NRVAL(r(b))
{
  NRCODE(r -= sl;)
  NONRCODE(Scl<T> r(b); r -= sl; return r;)
}


inline Scl<T>
operator * (const Scl<T>& b, long sl)  NRVAL(r(b))
{
  NRCODE(r *= sl;)
  NONRCODE(Scl<T> r(b); r *= sl; return r;)
}


inline Scl<T>
operator / (const Scl<T>& b, long sl)  NRVAL(r(b))
{
  NRCODE(r /= sl;)
  NONRCODE(Scl<T> r(b); r /= sl; return r;)
}


inline Scl<T>
operator % (const Scl<T>& b, long sl)  NRVAL(r(b))
{
  NRCODE(r %= sl;)
  NONRCODE(Scl<T> r(b); r %= sl; return r;)
}

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

inline Scl<T>
operator + (long sl, const Scl<T>& b)
{
  return b + sl;
}

inline Scl<T>
operator * (long sl, const Scl<T>& b)
{
  return b * sl;
}

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

inline Scl<T>
Div2Exp (const Scl<T>& b, unsigned long exponent_of_2) NRVAL(r(b))
{
  NRCODE(r.Div2Exp(exponent_of_2);)
  NONRCODE(Scl<T> r(b); r.Div2Exp(exponent_of_2); return r;)
}

inline Scl<T>
GCD (const Scl<T>& b1, const Scl<T>& b2) NRVAL(r(b1))
{
  NRCODE(r.GCD(b2);)
  NONRCODE(Scl<T> r(b1); r.GCD(b2); return r;)
}

inline Scl<T>
gcd (const Scl<T>& b1, const Scl<T>& b2) NRVAL(r(b1))
{
  NRCODE(r.GCD(b2);)
  NONRCODE(Scl<T> r(b1); r.GCD(b2); return r;)
}

inline Scl<T>
Mod2Exp (const Scl<T>& b, unsigned long exponent_of_2) NRVAL(r(b))
{
  NRCODE(r.Mod2Exp(exponent_of_2);)
  NONRCODE(Scl<T> r(b); r.Mod2Exp(exponent_of_2); return r;)
}

inline Scl<T>
Mul2Exp (const Scl<T>& b, unsigned long exponent_of_2) NRVAL(r(b))
{
  NRCODE(r.Mul2Exp(exponent_of_2);)
  NONRCODE(Scl<T> r(b); r.Mul2Exp(exponent_of_2); return r;)
}

inline Scl<T>
PowMod  (const Scl<T>& b, const Scl<T>& exp, const Scl<T>& m) NRVAL(r(b))
{
  NRCODE(r.PowMod(exp, m);)
  NONRCODE(Scl<T> r(b); r.PowMod(exp,m); return r;)
}

inline Scl<T>
PowMod  (const Scl<T>& b, unsigned long exp, const Scl<T>& m) NRVAL(r(b))
{
  NRCODE(r.PowMod(exp, m);)
  NONRCODE(Scl<T> r(b); r.PowMod(exp,m); return r;)
}

inline Scl<T>
abs (const Scl<T>& b) NRVAL(r(b))
{
  NRCODE(r.abs();)
  NONRCODE(Scl<T> r(b); r.abs(); return r;)  
}

// inline Scl<T>
// negate (const Scl<T>& b) NRVAL(r(b))
// {
//   NRCODE(r.negate();)
//   NONRCODE(Scl<T> r(b); r.negate(); return r;)  
// }

inline Scl<T>
pow (const Scl<T>& base, unsigned long exp) NRVAL(r(base))
{
  NRCODE(r.pow(exp);)
  NONRCODE(Scl<T> r(base); r.pow(exp); return r;)
}

inline Scl<T>
quo (const Scl<T>& dividend, const Scl<T>& divisor) NRVAL(r(dividend))
{
  NRCODE(r.quo(divisor);)
  NONRCODE(Scl<T> r(dividend); r.quo(divisor); return r;)
}

inline Scl<T>
quo (const Scl<T>& dividend, unsigned long divisor) NRVAL(r(dividend))
{
  NRCODE(r.quo(divisor);)
  NONRCODE(Scl<T> r(dividend); r.quo(divisor); return r;)
}

inline Scl<T>
rem (const Scl<T>& dividend, const Scl<T>& divisor) NRVAL(r(dividend))
{
  NRCODE(r.rem(divisor);)
  NONRCODE(Scl<T> r(dividend); r.rem(divisor); return r;)
}

inline Scl<T>
sqrt (const Scl<T>& b) NRVAL(r(b))
{
  NRCODE(r.sqrt();)
  NONRCODE(Scl<T> r(b); r.sqrt(); return r;)
}

inline unsigned long
rem (const Scl<T>& b, unsigned long divisor)
{
  Scl<T> r(b);
  return r.rem(divisor);
}

inline Scl<T>
BigIntFactorial (unsigned long n) NRVAL(r)
{
  NRCODE(r.factorial(n);)
  NONRCODE(Scl<T> r; r.factorial(n); return r;)    
}


inline Scl<T>
BigIntPow (unsigned long base, unsigned long exp) NRVAL(r)
{
  NRCODE(r.pow(base, exp);)
  NONRCODE(Scl<T> r; r.pow(base, exp); return r;)
}

#undef NRVAL
#undef NRCODE
#undef NONRCODE

__END_NAMESPACE_SYNAPS


#endif // SYNAPS_ARITHM__BC_H_INCLUDED_