GMP.hpp

Go to the documentation of this file.

#ifndef mmx_gmp_kernel_hpp
#define mmx_gmp_kernel_hpp
#include <realroot/texp_kernelof.hpp>
#include <realroot/scalar_integer.hpp>
#include <realroot/scalar_rational.hpp>
#include <realroot/scalar_floating.hpp>
#include <realroot/scalar_extended_integer.hpp>
#include <realroot/scalar_extended_rational.hpp>
#include <realroot/Interval.hpp>

namespace mmx {

struct GMP {
  typedef double      ieee;
  typedef scalar<MPZ> integer;
  typedef scalar<MPQ> rational;
  typedef scalar<extended<MPZ> > extended_integer;
  typedef scalar<extended<MPQ> > extended_rational;
  typedef scalar<MPF> floating;
  typedef Interval<rational> interval_rational;
  typedef Interval<extended_rational> interval_extended_rational;
};

  typedef GMP::rational QQ;
  typedef GMP::integer  ZZ;
  typedef GMP::floating RR;

  typedef GMP::extended_rational EQQ;
  typedef GMP::extended_integer  EZZ;

  // pdobrowolski: should these go to a separate texp_* file?
  template<class type> struct is_extended { typedef false_t T; };
  template<> struct is_extended<EQQ> { typedef true_t T; };
  template<> struct is_extended<EZZ> { typedef true_t T; };

  namespace texp {
    template<class R> struct kernelof_;
    template<> struct kernelof_<GMP::ieee>     { typedef GMP T; };
    template<> struct kernelof_<GMP::integer>  { typedef GMP T; };
    template<> struct kernelof_<GMP::rational> { typedef GMP T; };
    template<> struct kernelof_<GMP::floating> { typedef GMP T; };

    template<> struct kernelof_<GMP::extended_integer>  { typedef GMP T; };
    template<> struct kernelof_<GMP::extended_rational> { typedef GMP T; };
  }

  inline int  sign(const QQ& a) { return mpq_sgn(&a.rep()); }

  inline bool with_plus_sign(const ZZ & c) { return (c>0);}
  inline bool with_plus_sign(const QQ & c) { return (c>0);}
  inline bool with_plus_sign(const RR & c) { return (c>0);}
  
  inline double to_double(const RR& a)   { return mpf_get_d(&a.rep()); }
  inline double as_double(const RR& a)   { return mpf_get_d(&a.rep()); }

  inline double to_XT(const ZZ& a) { return as_double(a); }
  inline double to_XT(const QQ& a) { return as_double(a); }
  inline double to_XT(const RR& a) { return as_double(a); }

//----------------------------------------------------------------------
inline void bit_precision(unsigned long l) {
  mpf_set_default_prec(l);
}

//----------------------------------------------------------------------
inline void 
bit_precision( RR b, unsigned long l){
  mpf_set_prec(&b.rep(),l);
}

inline double gcd (const double&, const double&) { return 1; }


//======================================================================
namespace let
{
  inline void assign(ZZ& x, const QQ& r) { x= numerator(r)/denominator(r);}
  inline void assign(QQ& x, const ZZ& r) { mpq_set_z(&x.rep(), &r.rep());}
  inline void assign(RR& r, const ZZ& x) { mpf_set_z(&r.rep(),&x.rep());}

  inline void assign(mpz_t r, const ZZ& x)   { mpz_set(r, &x.rep());}
  inline void assign(mpf_t r, const ZZ& x)   { mpf_set_z(r,&x.rep());}
  inline void assign(mpf_t r, const RR& x)   { mpf_set(r,&x.rep());}
  inline void assign(mpf_t r, const QQ& x)
  {
    //    RR u=RR(numerator(x))/RR(denominator(x)); mpf_set(r,&u.rep());
    mpf_t n, d; mpf_init (n); mpf_init (d);
    assign (n, numerator (x));
    assign (d, denominator (x));
    mpf_div (r, n, d);
    mpf_clear (n);
    mpf_clear (d);
  }

}

template<class T, class U> struct as_helper;

template<> struct as_helper<QQ,ZZ>     { static inline QQ cv(const ZZ& x)     
  { QQ r;mpq_set_z(&r.rep(),&x.rep()); return r; } };

template<> struct as_helper<ZZ,QQ>     { static inline ZZ cv(const QQ& r) 
  {  return ( numerator(r)/denominator(r) ); } };

template<> struct as_helper<QQ,RR>     { static inline QQ cv(const RR& x)     
  { QQ r;mpq_set_f(&r.rep(),&x.rep()); return r; } };

template<> struct as_helper<RR,QQ>     { static inline RR cv(const QQ& x)     
  { RR r;mpf_set_q(&r.rep(),&x.rep()); return r; } };

template<> struct as_helper<ZZ,unsigned> { static inline ZZ cv(const unsigned& x)     
  { ZZ r;mpz_set_ui(&r.rep(),x); return r; } };

template<> struct as_helper<ZZ,double> { static inline ZZ cv(const double& x)     
  { ZZ r;mpz_set_d(&r.rep(),x); return r; } };

template<> struct as_helper<RR,double> { static inline RR cv(const double& x)     
  { RR r;mpf_set_d(&r.rep(),x); return r; } };

template<> struct as_helper<double,RR> { static inline double cv(const RR& x)     
  { return mpf_get_d(&x.rep()); } };

template<> struct as_helper<double,mpf_t> { static inline double cv(const mpf_t& x)     
  { return mpf_get_d(x); } };


template<> struct as_helper<Interval<double>,QQ> { 
  static inline Interval<double> cv(const QQ& x) { 
    double l,u;
    {
      numerics::rounding<double> rnd( numerics::rounding<double>::rnd_up() );
      u=as<double>(x);
    }
    {
      numerics::rounding<double> rnd( numerics::rounding<double>::rnd_dw() );
      l=as<double>(x);
    }

    return Interval<double>(l,u);
  } 
};

// template<> struct as_helper<double,ZZ> { static inline double cv(const ZZ& x) { return x.get_d(); } };
// template<> struct as_helper<double,RR> { static inline double cv(const RR& x) { return x.get_d(); } };


  inline QQ to_FT(const ZZ& a, const ZZ& b = 1) {
    QQ r,d; let::assign(r,a); let::assign(d,b);
    r /=d;
    //      r.canonicalize();
    return r;
  }
  inline QQ to_FT(const QQ& a, const QQ& b = 1) { return a/b; }
//====================================================================
} //namespace mmx
#endif