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

Go to the documentation of this file.

#ifndef realroot_ARITHM_GMP_H
#define realroot_ARITHM_GMP_H

#include <realroot/config.hpp>
#include <realroot/rounding_mode.hpp>
#include <realroot/Interval.hpp>
#include <realroot/sign.hpp>

#include <gmpxx.h>


#ifdef HAVE_MPFR
#include <mpfr.h>
#endif
namespace mmx {

struct GMP
{
  typedef mpz_class integer;
  typedef mpq_class rational;
  typedef mpf_class floating;
  typedef double    ieee;
  
};


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


inline long int bit_size(const ZZ& z)  {return mpz_sizeinbase(z.get_mpz_t(),2);}

namespace meta {

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

  template<class C> struct hasgcd_;
  template<> struct hasgcd_<ZZ> { typedef meta::true_t T; };
  
  template<class C> struct structureof_;
  template<> struct structureof_< GMP::rational > { typedef structure::scalar T; };
  template<> struct structureof_< GMP::integer  > { typedef structure::scalar T; };
  template<> struct structureof_< GMP::floating > { typedef structure::scalar T; };
}
  
  inline ZZ gcd(const ZZ & a, const ZZ & b)
    {
      mpz_class r;
      mpz_gcd(r.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
      return r;
    }
  
  inline ZZ lcm(const ZZ & a, const ZZ & b)  { return (a*b)/gcd(a,b); }


  inline int sign(const QQ& a) { return sgn(a); }
  inline int sign(const RR& a) { return sgn(a); }

  inline QQ Size(const QQ & r) { return abs(r); };
  inline ZZ Size( const ZZ & z ) { return abs(z); };

  inline ZZ size( const ZZ & z ) { return abs(z); };
  inline QQ size(const QQ & r) { return abs(r); };
  
  
  inline int compare(const QQ& a, const QQ& b) { return cmp(a, b); }
  inline int compare(const RR& a, const RR& b) { return cmp(a, b); }

  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 QQ& a) { return a.get_d(); }
  inline double to_double( const ZZ & z ) { return z.get_d(); };
  inline double to_double(const RR& a) { return a.get_d(); }

  inline QQ to_FT(const ZZ& a, const ZZ& b = 1)
    {
      QQ r(a, b);
      r.canonicalize();
      return r;
    }
  inline QQ to_FT(const QQ& a, const QQ& b = 1) { return a; }

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

  inline ZZ numerator(const QQ& a)   { return a.get_num(); }
  inline ZZ denominator(const QQ& a) { return a.get_den(); }

  inline ZZ pow(const ZZ& a, unsigned  n)
    {
      mpz_class r;
      mpz_pow_ui(r.get_mpz_t(), a.get_mpz_t(), n);
      return r;
    }
  inline ZZ isqrt(const ZZ& a)    { return sqrt(a) + 1; }
  inline RR fracpart(const RR& r) { return r - trunc(r); }


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

//----------------------------------------------------------------------
inline void Precision( RR b, unsigned long l)
{
    mpf_set_prec(b.get_mpf_t(),l);

    }

/*
template<class X>
inline void precision( unsigned nb )
{
  typedef kernelof<X>::T K;
  K::precision(nb);
};
*/
//======================================================================

namespace let
{
  inline void assign(QQ & x, char *s) { mpq_set_str(x.get_mpq_t(), s, 10);}

  inline void assign(ZZ& z, char * s)   { mpz_set_str(z.get_mpz_t(), s, 10); }
  inline void assign(ZZ& z, int n)      { mpz_set_si(z.get_mpz_t(),n); }
  inline void assign(ZZ& z, double d )  { z = (int)d; };
  inline void assign(ZZ& x, const ZZ& r)      { mpz_set(x.get_mpz_t(),r.get_mpz_t()); }
  inline void assign(int& x, const ZZ& r)       { x = mpz_get_si(r.get_mpz_t()); }
  inline void assign(long int& x, const ZZ& r)  { x = mpz_get_si(r.get_mpz_t()); }
  inline void assign( double & r, const ZZ & z ) { r = z.get_d(); };
  
  inline void assign(ZZ& x, const QQ& r) { x= r.get_num(); x/=r.get_den();}
  inline void assign(QQ& x, const ZZ& r) { mpq_set_z(x.get_mpq_t(), r.get_mpz_t());}
  inline void assign(QQ & q, double d ) { q = d; };
  inline void assign(RR & r, double d ) { r = d; };
  inline void assign(double & d, const QQ & x ) 
  {
#ifdef HAVE_MPFR 
    mp_rnd_t round=GMP_RNDU;
    //    if(numerics::get_rnd()==numerics::rnd_up())
    //      round=GMP_RNDU;

    if(numerics::get_rnd()==numerics::rnd_dw())
      round=GMP_RNDD;

    mpfr_t r; mpfr_init2(r,53);
    mpfr_set_q(r, x.get_mpq_t(), round);
    d = mpfr_get_d(r,round);
#else
       d = x.get_d(); 
#endif
  }

  inline void assign(double & r, const RR & z ) { r = z.get_d(); };  
  inline void assign(RR& r, char* s)        { mpf_set_str(r.get_mpf_t(),s,10); }
  inline void assign(RR& x, const RR& r)    { mpf_set(x.get_mpf_t(),r.get_mpf_t()); }

  inline double convert( const QQ & q, meta::As<double> ) { return q.get_d(); };
  inline double convert( const ZZ & z, meta::As<double> ) { return z.get_d(); };
  inline double convert( const RR & r, meta::As<double> ) { return r.get_d(); };

  inline void assign( Interval<double>& a, 
                      const QQ & b )
    {
      double  mn,mx;
      {
        numerics::rounding<double> rnd( numerics::rnd_up() );
        assign(mx,b);
      }
      {
        numerics::rounding<double> rnd( numerics::rnd_dw() );
        assign(mn,b);
      }
      
      if( mx < mn ) std::swap(mn,mx);
      QQ MN(mn);
      QQ MX(mx);
      new (&a) Interval<double>(mn,mx);
    };
  inline void assign( Interval<double>& a, 
                      const QQ & u,
                      const QQ & v )
    {
      double  mn,mx;
      {
        numerics::rounding<double> rnd( numerics::rnd_up() );
        assign(mx,v);
      }
      {
        numerics::rounding<double> rnd( numerics::rnd_dw() );
        assign(mn,u);
      }
      new (&a) Interval<double>(mn,mx);
    };


  inline void assign( Interval<QQ>& a, 
                      const QQ & u,
                      const QQ & v )
  {
    a.lower()=u;
    a.upper()=v;
  }
}

template<class T, class U> struct cast_helper;

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

template<> struct cast_helper<double,ZZ> { static inline double cv(const ZZ& x) { return x.get_d(); } };
template<> struct cast_helper<double,QQ> { static inline double cv(const QQ& x) { double r; let::assign(r,x); return r; } };
template<> struct cast_helper<double,RR> { static inline double cv(const RR& x) { return x.get_d(); } };

//====================================================================
} //namespace mmx
#endif