/Users/mourrain/Devel/mmx/realroot/include/realroot/scalar_integer.hpp

Go to the documentation of this file.

/*********************************************************************
*      This file is part of the source code of ALP software.         *
*      (C) B. Mourrain, GALAAD, INRIA                                *
**********************************************************************
History: 
$Id: scalar_mpz.H,v 1.1.1.1 2006/10/06 08:01:40 trebuche Exp $
**********************************************************************/
#ifndef _scalar_integer_hpp
#define _scalar_integer_hpp
//----------------------------------------------------------------------
#include <iostream>
#include <string>
#include <gmp.h>
#include <realroot/scalar.hpp>
//======================================================================
namespace mmx {
//----------------------------------------------------------------------
typedef MP_INT MPZ;
//----------------------------------------------------------------------
//  template<> inline
//  shared_object<MPZ>::rep::~rep()           
//  {
//    //COUT <<"Delete mpz"<<endl;
//    mpz_clear(&obj);
//  }
//----------------------------------------------------------------------
template<> inline scalar<MPZ>::~scalar ( )      
{ 
  //COUT <<"Delete mpz"<<endl;
  mpz_clear(&rep());
}
//----------------------------------------------------------------------
template<> inline
void scalar<MPZ>::init ( )                   
{ 
  mpz_init(&rep());
}
//------------------- scalar<T> (...) -------------------------------------
template<> inline
scalar<MPZ>::scalar () //: data()            
{
  mpz_init(&rep());
}
//----------------------------------------------------------------------
template<> inline
scalar<MPZ>::scalar (signed long int sl): data()  {mpz_init_set_si(&rep(),sl);}
//----------------------------------------------------------------------
template<> inline
scalar<MPZ>::scalar (unsigned long int ul): data() {mpz_init_set_ui(&rep(),ul);}
//----------------------------------------------------------------------
template<> inline
scalar<MPZ>::scalar (int si): data()               {mpz_init_set_si(&rep(), si);}
//----------------------------------------------------------------------
template<> inline
scalar<MPZ>::scalar (const char *s, unsigned int base): data()
{
  if (mpz_init_set_str(&rep(), s, base))
    std::cerr << "scalar<MPZ>: The string " << s 
              << " is not a valid number in base " << base << std::endl;
}
//----------------------------------------------------------------------
template<> inline
scalar<MPZ>::scalar (const scalar<MPZ> & rhs)
{
  //COUT <<"Init mpz"<<endl;
  //COUNT<MPZ>('c');
  mpz_init_set(&rep(), &rhs.rep());      
}
//{ data=b.data;}
//----------------------------------------------------------------------
template<> inline 
scalar<MPZ>& scalar<MPZ>::operator =  (const scalar<MPZ>& rhs)
{
  //COUNT<MPZ>('=');
  if (this != &rhs) mpz_set(&rep(), &rhs.rep()); 
  return *this;
}
template<> inline 
scalar<MPZ>& scalar<MPZ>::operator =  (unsigned rhs)
{
  //COUNT<MPZ>('=');
  mpz_set_ui(&rep(), rhs); 
  return *this;
}
// == ------------------------------------------------------------------
template<> inline 
bool scalar<MPZ>::operator == (const scalar<MPZ>& rhs) const
{
  return mpz_cmp(&rep(), &rhs.rep()) == 0;
}

template<> inline 
bool scalar<MPZ>::operator == (long sl) const {
  return mpz_cmp_si(&rep(), sl) == 0;
}

template<> inline 
bool scalar<MPZ>::operator == (int si) const {
  return mpz_cmp_si(&rep(), (long) si) == 0;
}

template<> inline 
bool scalar<MPZ>::operator == (unsigned long ul) const
{
  return mpz_cmp_ui(&rep(), ul) == 0;
}

// != ------------------------------------------------------------------
template<> inline 
bool scalar<MPZ>::operator != (const scalar<MPZ>& rhs) const
{
  return mpz_cmp(&rep(), &rhs.rep()) != 0;
}

template<> inline 
bool scalar<MPZ>::operator != (long sl) const {
  return mpz_cmp_si(&rep(), sl) != 0;
}

template<> inline
bool scalar<MPZ>::operator != (int si) const {
  return mpz_cmp_si(&rep(), (long) si) != 0;
}

template<> inline 
bool scalar<MPZ>::operator != (unsigned long ul) const
{
  return mpz_cmp_ui(&rep(), ul) != 0;
}
// > -------------------------------------------------------------------
template<> inline 
bool scalar<MPZ>::operator > (const scalar<MPZ>& rhs) const
{
  return mpz_cmp(&rep(), &rhs.rep()) > 0;
}

template<> inline 
bool scalar<MPZ>::operator > (long sl) const {
  return mpz_cmp_si(&rep(), sl) > 0;
}

template<> inline 
bool scalar<MPZ>::operator > (int si) const {
  return mpz_cmp_si(&rep(), (long) si) > 0;
}

template<> inline 
bool scalar<MPZ>::operator > (unsigned long ul) const
{
  return mpz_cmp_ui(&rep(), ul) > 0;
}
// <= --------------------------------------------------------------------
template<> inline 
bool scalar<MPZ>::operator >= (const scalar<MPZ>& rhs) const
{
  return mpz_cmp(&rep(), &rhs.rep()) >= 0;
}

template<> inline 
bool scalar<MPZ>::operator >= (long sl) const {
  return mpz_cmp_si(&rep(), sl) >= 0;
}

template<> inline 
bool scalar<MPZ>::operator >= (int si) const
{
  return mpz_cmp_si(&rep(), (long) si) >= 0;
}

template<> inline 
bool scalar<MPZ>::operator >= (unsigned long ul) const
{
  return mpz_cmp_ui(&rep(), ul) >= 0;
}

// < -------------------------------------------------------------------
template<> inline 
bool scalar<MPZ>::operator < (const scalar<MPZ>& rhs) const
{
  return mpz_cmp(&rep(), &rhs.rep()) < 0;
}

template<> inline 
bool scalar<MPZ>::operator < (long sl) const {
  return mpz_cmp_si(&rep(), sl) < 0;
}

template<> inline 
bool scalar<MPZ>::operator < (int si) const
{
  return mpz_cmp_si(&rep(), (long) si) < 0;
}

template<> inline 
bool scalar<MPZ>::operator < (unsigned long ul) const
{
  return mpz_cmp_ui(&rep(), ul) < 0;
}

// <= ------------------------------------------------------------------

template<> inline 
bool scalar<MPZ>::operator <= (const scalar<MPZ>& rhs) const
{
  return mpz_cmp(&rep(), &rhs.rep()) <= 0;
}

template<> inline 
bool scalar<MPZ>::operator <= (long sl) const {
  return mpz_cmp_si(&rep(), sl) <= 0;
}

template<> inline 
bool scalar<MPZ>::operator <= (int si) const
{
  return mpz_cmp_si(&rep(), (long) si) <= 0;
}

template<> inline 
bool scalar<MPZ>::operator <= (unsigned long ul) const
{
  return mpz_cmp_ui(&rep(), ul) <= 0;
}

// = -------------------------------------------------------------------
template<> inline 
scalar<MPZ>&  scalar<MPZ>::operator =  (unsigned long ul)
{
  mpz_set_ui(&rep(), ul); return *this;
}

template<> inline 
scalar<MPZ>& scalar<MPZ>::operator =  (long sl)
{
  mpz_set_si(&rep(), sl);  return *this;
}

template<> inline 
scalar<MPZ>& scalar<MPZ>::operator = (int ul)
{
  mpz_init_set_si(&(this->rep()), ul);  return *this;
}
// += -------------------------------------------------------------------
template<> inline 
scalar<MPZ>& scalar<MPZ>::operator += (const scalar<MPZ>& rhs)
{
  mpz_add(&rep(), &rep(), &rhs.rep());  return *this;
}

template<> inline 
scalar<MPZ>& scalar<MPZ>::operator += (unsigned long ul)
{
  mpz_add_ui(&rep(), &rep(), ul); return *this;
}

template<> inline 
scalar<MPZ>& scalar<MPZ>::operator += (long sl)
{
  if (sl >= 0)
    mpz_add_ui(&rep(), &rep(), (unsigned long) sl);
  else
    mpz_sub_ui(&rep(), &rep(), ((unsigned long) -sl));
  return *this;
}

template<> inline 
scalar<MPZ>& scalar<MPZ>::operator += (int ul)
{
  *this += scalar<MPZ>(ul);  return *this;
}

// -= -------------------------------------------------------------------
template<> inline 
scalar<MPZ>& scalar<MPZ>::operator -= (const scalar<MPZ>& rhs)
{
  mpz_sub(&rep(), &rep(), &rhs.rep()); return *this;
}

template<> inline 
scalar<MPZ>& scalar<MPZ>::operator -= (unsigned long ul)
{
  mpz_sub_ui(&rep(), &rep(), ul); return *this;
}

template<> inline 
scalar<MPZ>& scalar<MPZ>::operator -= (long sl)
{
  if (sl >= 0)
    mpz_sub_ui(&rep(), &rep(), (unsigned long) sl);
  else
    mpz_add_ui(&rep(), &rep(), (unsigned long) sl);
  return *this;
}

template<> inline 
scalar<MPZ>& scalar<MPZ>::operator -= (int ul)
{
  *this -= scalar<MPZ>(ul); return *this;
}
// *= -----------------------------------------------------------------
template<> inline 
scalar<MPZ>& scalar<MPZ>::operator *= (const scalar<MPZ>& rhs)
{
  mpz_mul(&rep(), &rep(), &rhs.rep()); return *this;
}
template<> inline 
scalar<MPZ>& scalar<MPZ>::operator *= (unsigned long ul)
{
  mpz_mul_ui(&rep(), &rep(), ul); return *this;
}
template<> inline 
scalar<MPZ>& scalar<MPZ>::operator *= (long sl)
{
  if (sl >= 0)
    mpz_mul_ui(&rep(), &rep(), (unsigned long) sl);
  else
    {
      rep()._mp_size = -rep()._mp_size;
      mpz_mul_ui(&rep(), &rep(), ((unsigned long) -sl));
    }
  return *this; 
}
template<> inline 
scalar<MPZ>& scalar<MPZ>::operator *= (int ul)
{
   if (ul >= 0)
     mpz_mul_ui(&rep(), &rep(), (unsigned long) ul);
   else {
     mpz_mul_ui(&rep(), &rep(), (unsigned long) (-ul));
     mpz_neg(&rep(), &rep());
   }
   return *this;
}
// /= ------------------------------------------------------------------
template<> inline 
scalar<MPZ>& scalar<MPZ>::operator /= (const scalar<MPZ>& rhs)
{
  mpz_div(&rep(), &rep(), &rhs.rep());
  //mpz_divexact(&rep(), &rep(), &rhs.rep()); 
  return *this;
}
template<> inline 
scalar<MPZ>& scalar<MPZ>::operator /= (unsigned long ul)
{
  mpz_div_ui(&rep(), &rep(), ul); return *this;
}

template<> inline 
scalar<MPZ>& scalar<MPZ>::operator /= (long sl)
{
  if (sl >= 0)
    mpz_div_ui(&rep(), &rep(), (unsigned long) sl);
  else
    {
      rep()._mp_size = -rep()._mp_size;
      mpz_div_ui(&rep(), &rep(), ((unsigned long) -sl));
    }
  return *this;
}

// %= ------------------------------------------------------------------
template<> inline 
scalar<MPZ>& scalar<MPZ>::operator %= (const scalar<MPZ>& rhs)
{
  mpz_mod(&rep(), &rep(), &rhs.rep()); return *this;
}

template<> inline 
scalar<MPZ>& scalar<MPZ>::operator %= (unsigned long ul)
{
  mpz_mod_ui(&rep(), &rep(), ul); return *this;
}

//======================================================================
// ARITHMETIC OPERATORS
//======================================================================
inline
scalar<MPZ> operator +(const scalar<MPZ>& a1, const scalar<MPZ>& a2) 
{
  scalar<MPZ> result;
  mpz_add(&result.rep(), &a1.rep(), &a2.rep());
  return result;
}
//----------------------------------------------------------------------
inline
scalar<MPZ> operator -(const scalar<MPZ>& a1, const scalar<MPZ>& a2) 
{
   scalar<MPZ> result;
   mpz_sub(&result.rep(), &a1.rep(), &a2.rep());
   return result;
 }

inline
scalar<MPZ> operator -(const scalar<MPZ>& a1) 
{
  scalar<MPZ> r;  mpz_neg(&r.rep(), &a1.rep()); return r;
}
//----------------------------------------------------------------------
inline
scalar<MPZ> operator *(const scalar<MPZ>& a1, const scalar<MPZ>& a2) 
{
  scalar<MPZ> result;
  mpz_mul(&result.rep(), &a1.rep(), &a2.rep());
  return result;
}

//--------------------------------------------------------------------
inline 
scalar<MPZ> operator /(const scalar<MPZ>& a1, const scalar<MPZ>& a2) 
{
  scalar<MPZ> result;
  mpz_div(&result.rep(), &a1.rep(), &a2.rep());
  return result;
}
//======================================================================
template<> inline 
void scalar<MPZ>::operator ++ ( )
{
  mpz_add_ui(&rep(), &rep(), 1);
}
//----------------------------------------------------------------------
template<> inline 
void scalar<MPZ>::operator -- ( )
{
  mpz_sub_ui(&rep(), &rep(), 1);
}
//======================================================================
inline void convert(scalar<MPZ>& n, char *s) 
{
  mpz_init_set_str(&n.rep(), s, 10);
}
//======================================================================
// INPUT OUTPUT
//======================================================================
inline
std::ostream& operator << (std::ostream& os, const scalar<MPZ>& b)
{
  mpz_out_str(stdout, 10, &b.rep());
  return os;
}
//----------------------------------------------------------------------
inline
std::istream& operator >> (std::istream& is, scalar<MPZ>& b)
{
  std::string s;is >> s;
  mpz_init_set_str(&b.rep(),s.c_str(), 10);
  return is;
}
//======================================================================
// SPECIAL FUNCTIONS
//======================================================================
template<> inline void
scalar<MPZ>::Div2Exp (unsigned long exponent_of_2)
{
  mpz_div_2exp(&rep(), &rep(), exponent_of_2);
}

template<> inline 
void scalar<MPZ>::GCD (const scalar<MPZ>& b2)
{
    mpz_gcd (&rep(), &rep(), &b2.rep());
}


template<> inline void
scalar<MPZ>::Mod2Exp (unsigned long exponent_of_2)
{
  mpz_mod_2exp(&rep(), &rep(), exponent_of_2);
}

template<> inline void
scalar<MPZ>::Mul2Exp (unsigned long exponent_of_2)
{
  mpz_mul_2exp(&rep(), &rep(), exponent_of_2);
}

template<> inline void
scalar<MPZ>::PowMod (const scalar<MPZ>& exp, const scalar<MPZ>& m)
{
  mpz_powm(&rep(), &rep(), &exp.rep(), &m.rep());
}

template<> inline void
scalar<MPZ>::PowMod (unsigned long exp, const scalar<MPZ>& m)
{
  mpz_powm_ui(&rep(), &rep(), exp, &m.rep());
}

template<> inline void
scalar<MPZ>::abs ( )
{
  if (rep()._mp_size < 0)
    rep()._mp_size = - rep()._mp_size;
}

template<> inline void
scalar<MPZ>::factorial (unsigned long n)
{
  mpz_fac_ui(&rep(), n);
}

template<> inline 
scalar<MPZ>& scalar<MPZ>::negate ( )
{
  rep()._mp_size = - rep()._mp_size; return *this;
}

template<> inline 
scalar<MPZ>& scalar<MPZ>::pow (unsigned long exp)
{
  mpz_pow_ui(&rep(), &rep(), exp); return *this;
}

template<> inline 
scalar<MPZ>& scalar<MPZ>::pow (unsigned long base, unsigned long exp)
{
  mpz_ui_pow_ui(&rep(), base, exp); return *this;
}

template<> inline void
scalar<MPZ>::quo (const scalar<MPZ>& divisor)
{
  mpz_mdiv (&rep(), &rep(), &divisor.rep());
}

template<> inline void
scalar<MPZ>::quo (unsigned long divisor)
{
  mpz_mdiv_ui(&rep(), &rep(), divisor);
}

template<> inline void
scalar<MPZ>::rem (const scalar<MPZ>& divisor)
{
  mpz_mmod(&rep(), &rep(), &divisor.rep());
}

template<> inline unsigned long
scalar<MPZ>::rem (unsigned long divisor)
{
  return mpz_mmod_ui(&rep(), &rep(), divisor);
}

template<> inline void
scalar<MPZ>::sqrt ( )
{
  mpz_sqrt(&rep(), &rep());
}

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

inline unsigned long
QuoRem (scalar<MPZ>& q, scalar<MPZ>& r, const scalar<MPZ>& divdend, unsigned long divisor)
{
  return mpz_mdivmod_ui(&q.rep(), &r.rep(), &divdend.rep(), divisor);
}

inline void
DivMod (scalar<MPZ>& q, scalar<MPZ>& r, const scalar<MPZ>& divdend, const scalar<MPZ>& divisor)
{
  mpz_divmod(&q.rep(), &r.rep(), &divdend.rep(), &divisor.rep());
}

inline void
DivMod (scalar<MPZ>& q, scalar<MPZ>& r, const scalar<MPZ>& divdend, unsigned long divisor)
{
  mpz_divmod_ui(&q.rep(), &r.rep(), &divdend.rep(), divisor);
}

inline 
void ExtGCD (scalar<MPZ>& gcd, scalar<MPZ>& a, scalar<MPZ>& b, 
             const scalar<MPZ>& x, const scalar<MPZ>& y)
{
  mpz_gcdext(&gcd.rep(), &a.rep(), &b.rep(), &x.rep(), &y.rep());
}

inline 
void HalfExtGCD (scalar<MPZ>& gcd, scalar<MPZ>& a, 
                 const scalar<MPZ>& x, const scalar<MPZ>& y)
{
    mpz_gcdext(&gcd.rep(), &a.rep(), 0, &x.rep(), &y.rep());
}

inline unsigned long BigIntToUL (const scalar<MPZ>& b)
{
  return mpz_get_ui(&b.rep());
}
//-------------------------- BigIntToSL ---------------------------------
inline signed long BigIntToSL (const scalar<MPZ>& b)
{
  return mpz_get_si(&b.rep());
}
//-------------------------- log(...) -----------------------------------
inline size_t log (const scalar<MPZ>& b)
{
  return mpz_size(&b.rep());
}

inline size_t log (const scalar<MPZ>& b, int base)
{
  assert(base >= 2 && base <= 36);
  return mpz_sizeinbase(&b.rep(), base);
}
//-------------------------- sign ---------------------------------------------
inline int sign (const scalar<MPZ>& b)
{
  if (b.rep()._mp_size == 0)
    return 0;
  else
    return b.rep()._mp_size > 0 ? 1 : -1;
}
//-------------------------- compare (...) ------------------------------------
inline int compare (const scalar<MPZ>& b1, const scalar<MPZ>& b2)
{
  return mpz_cmp(&b1.rep(), &b2.rep());
}

inline int compare (const scalar<MPZ>& b, unsigned long ul)
{
  return mpz_cmp_ui(&b.rep(), ul);
}

inline int compare (const scalar<MPZ>& b, long sl)
{
  return mpz_cmp_si(&b.rep(), sl);
}
//-------------------------- IsPositive etc. ----------------------------------
inline bool IsPositive (const scalar<MPZ>& b) {
  return b.rep()._mp_size > 0;
}

inline bool IsNegative (const scalar<MPZ>& b) {
  return b.rep()._mp_size < 0;
}

inline bool IsZero (const scalar<MPZ>& b) {
  return b.rep()._mp_size == 0;
}

inline bool IsOne(const scalar<MPZ>& b) {
  return b.rep()._mp_size == 1 && b.rep()._mp_d[0] == 1;
}

inline bool IsMinusOne (const scalar<MPZ>& b) {
  return b.rep()._mp_size == -1 && b.rep()._mp_d[0] == 1;
}

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

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

inline bool IsPerfectSquare (const scalar<MPZ>& b)
{
  return mpz_perfect_square_p(&b.rep());
}

inline bool IsProbablyPrime (const scalar<MPZ>& b, int reps /* = 25 */)
{
  return mpz_probab_prime_p(&b.rep(), reps);
}
//----------------------------------------------------------------------
inline
scalar<MPZ> operator << (const scalar<MPZ>& x, long int s) {
  scalar<MPZ> r; 
  if (s >= 0) mpz_mul_2exp (&r.rep(), &x.rep(), s);
  else mpz_div_2exp (&r.rep(), &x.rep(), -s);
  return r; 
}
//----------------------------------------------------------------------
inline
scalar<MPZ>& operator <<= (scalar<MPZ>& x, long int s) 
{
  if (s >= 0) mpz_mul_2exp (&x.rep(), &x.rep(), s);
  else mpz_div_2exp (&x.rep(), &x.rep(), -s);
  return x; 
}

inline scalar<MPZ> Size( const scalar<MPZ> & z ) { return abs(z); };

inline long int bit_size(const scalar<MPZ>& z)  {return mpz_sizeinbase(&z.rep(),2);}

inline scalar<MPZ> gcd(const scalar<MPZ> & a, const scalar<MPZ> & b)
{
  scalar<MPZ> r;
  mpz_gcd(&r.rep(), &a.rep(), &b.rep());
  return r;
}

inline scalar<MPZ> lcm(const scalar<MPZ> & a, const scalar<MPZ> & b)  { return (a*b)/gcd(a,b); }
//--------------------------------------------------------------------
inline scalar<MPZ> pow(const scalar<MPZ>& a, unsigned  n)
{
  scalar<MPZ> r;
  mpz_pow_ui(&r.rep(), &a.rep(), n);
  return r;
}

inline scalar<MPZ> isqrt(const scalar<MPZ>& a)    { return sqrt(a) + 1; }

//--------------------------------------------------------------------
namespace let
{
  inline void assign(scalar<MPZ>& z, char * s)   { mpz_set_str(&z.rep(), s, 10); }
  inline void assign(scalar<MPZ>& z, int n)      { mpz_set_si(&z.rep(),n); }
  inline void assign(scalar<MPZ>& z, double d )  { z = (int)d; };
  inline void assign(scalar<MPZ>& x, const scalar<MPZ>& r)      { mpz_set(&x.rep(),&r.rep()); }
  inline void assign(int& x, const scalar<MPZ>& r)       { x = mpz_get_si(&r.rep()); }
  inline void assign(long int& x, const scalar<MPZ>& r)  { x = mpz_get_si(&r.rep()); }
  inline void assign( double & r, const scalar<MPZ> & z ) { r = mpz_get_d(&z.rep()); }
  
}

//   inline double convert( const scalar<MPZ> & z, meta::As<double> ) { return z.get_d(); };
  inline double to_double(const scalar<MPZ>& z)   { return mpz_get_d(&z.rep()); }

  template<typename T,typename F> struct as_helper;
  template<> struct as_helper<double,scalar<MPZ> > { 
  static inline double cv(const scalar<MPZ>& x) {return mpz_get_d(&x.rep());}
};
} // namespace mmx
//======================================================================
#endif // //SCL_MPZ_H