synaps/arithm/gmp.h

#ifndef SYNAPS_ARITHM_GMP_H
#define SYNAPS_ARITHM_GMP_H

#include <synaps/init.h>
#include <synaps/base/type.h>
#include <synaps/arithm/let.h>
#include <synaps/arithm/traits.h>
#include <synaps/arithm/functions.h>
#include <synaps/arithm/rounding.h>
#include <synaps/arithm/Interval.h>
#include <synaps/arithm/traits/structureof.h>
#include <gmpxx.h>

__BEGIN_NAMESPACE_SYNAPS

struct GMPSet
{
  typedef mpz_class integer;
  typedef mpq_class rational;
  typedef mpf_class floating;
  
};

typedef kernel<GMPSet> GMP;

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


namespace arithm {

  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<> struct hasgcd_<ZZ> { typedef type::true_t T; };
  
  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(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 & r, const QQ & z ) { r = z.get_d(); };
  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, type::As<double> ) { return q.get_d(); };
  inline double convert( const ZZ & z, type::As<double> ) { return z.get_d(); };
  inline double convert( const RR & r, type::As<double> ) { return r.get_d(); };
  //  template<class K>
  inline void assign( Interval<double, 3>& 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,3>(mn,mx);
    };

}

__END_NAMESPACE_SYNAPS

#endif