synaps/arithm/gmp.h

00001 #ifndef SYNAPS_ARITHM_GMP_H
00002 #define SYNAPS_ARITHM_GMP_H
00003 
00004 #include <synaps/init.h>
00005 #include <synaps/base/type.h>
00006 #include <synaps/arithm/let.h>
00007 #include <synaps/arithm/traits.h>
00008 #include <synaps/arithm/functions.h>
00009 #include <synaps/arithm/rounding.h>
00010 #include <synaps/arithm/Interval.h>
00011 #include <synaps/arithm/traits/structureof.h>
00012 #include <gmpxx.h>
00013 
00014 __BEGIN_NAMESPACE_SYNAPS
00015 
00016 struct GMPSet
00017 {
00018   typedef mpz_class integer;
00019   typedef mpq_class rational;
00020   typedef mpf_class floating;
00021   
00022 };
00023 
00024 typedef kernel<GMPSet> GMP;
00025 
00026 typedef GMP::rational QQ;
00027 typedef GMP::integer  ZZ;
00028 typedef GMP::floating RR;
00029 
00030 
00031 namespace arithm {
00032 
00033   template<> struct kernelof_< GMP::rational > { typedef GMP T; };
00034   template<> struct kernelof_< GMP::integer  > { typedef GMP T; };
00035   template<> struct kernelof_< GMP::floating > { typedef GMP T; };
00036   template<> struct hasgcd_<ZZ> { typedef type::true_t T; };
00037   
00038   template<> struct structureof_< GMP::rational > { typedef structure::scalar T; };
00039   template<> struct structureof_< GMP::integer  > { typedef structure::scalar T; };
00040   template<> struct structureof_< GMP::floating > { typedef structure::scalar T; };
00041 }
00042 
00043   inline ZZ gcd(const ZZ & a, const ZZ & b)
00044     {
00045       mpz_class r;
00046       mpz_gcd(r.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
00047       return r;
00048     }
00049   
00050   inline ZZ lcm(const ZZ & a, const ZZ & b)  { return (a*b)/gcd(a,b); }
00051 
00052 
00053   inline int sign(const QQ& a) { return sgn(a); }
00054   inline int sign(const RR& a) { return sgn(a); }
00055 
00056   inline QQ Size(const QQ & r) { return abs(r); };
00057   inline ZZ Size( const ZZ & z ) { return abs(z); };
00058 
00059   inline ZZ size( const ZZ & z ) { return abs(z); };
00060   inline QQ size(const QQ & r) { return abs(r); };
00061   
00062   
00063   inline int compare(const QQ& a, const QQ& b) { return cmp(a, b); }
00064   inline int compare(const RR& a, const RR& b) { return cmp(a, b); }
00065 
00066   inline bool with_plus_sign(const ZZ & c) { return (c>0);}
00067   inline bool with_plus_sign(const QQ & c) { return (c>0);}
00068   inline bool with_plus_sign(const RR & c) { return (c>0);}
00069   
00070 
00071   inline double to_double(const QQ& a) { return a.get_d(); }
00072   inline double to_double( const ZZ & z ) { return z.get_d(); };
00073   inline double to_double(const RR& a) { return a.get_d(); }
00074 
00075   inline QQ to_FT(const ZZ& a, const ZZ& b = 1)
00076     {
00077       QQ r(a, b);
00078       r.canonicalize();
00079       return r;
00080     }
00081   inline QQ to_FT(const QQ& a, const QQ& b = 1) { return a; }
00082 
00083   inline double to_XT(const ZZ& a) {  return to_double(a);  }
00084   inline double to_XT(const QQ& a) { return to_double(a); }
00085   inline double to_XT(const RR& a) { return to_double(a); }
00086 
00087   inline ZZ numerator(const QQ& a) { return a.get_num(); }
00088   inline ZZ denominator(const QQ& a) { return a.get_den(); }
00089 
00090   inline ZZ pow(const ZZ& a, unsigned  n)
00091     {
00092       mpz_class r;
00093       mpz_pow_ui(r.get_mpz_t(), a.get_mpz_t(), n);
00094       return r;
00095     }
00096   inline ZZ isqrt(const ZZ& a) { return sqrt(a) + 1; }
00097   inline RR fracpart( const RR& r)  { return r - trunc(r); }
00098 
00099 
00100 
00101 //----------------------------------------------------------------------
00105 inline void Precision(unsigned long l)
00106 {
00107     mpf_set_default_prec(l);
00108 }
00109 
00110 //----------------------------------------------------------------------
00114 inline void Precision( RR b, unsigned long l)
00115 {
00116     mpf_set_prec(b.get_mpf_t(),l);
00117 
00118     }
00119 
00120 /*
00121 template<class X>
00122 inline void precision( unsigned nb )
00123 {
00124   typedef kernelof<X>::T K;
00125   K::precision(nb);
00126 };
00127 */
00128 //======================================================================
00129 
00130 namespace let
00131 {
00132   inline void assign(QQ & x, char *s) { mpq_set_str(x.get_mpq_t(), s, 10);}
00133 
00134   inline void assign(ZZ& z, char * s)   { mpz_set_str(z.get_mpz_t(), s, 10); }
00135   inline void assign(ZZ& z, int n)      { mpz_set_si(z.get_mpz_t(),n); }
00136   inline void assign(ZZ& z, double d )  { z = (int)d; };
00137   inline void assign(ZZ& x, const ZZ& r)      { mpz_set(x.get_mpz_t(),r.get_mpz_t()); }
00138   inline void assign(int& x, const ZZ& r)       { x = mpz_get_si(r.get_mpz_t()); }
00139   inline void assign(long int& x, const ZZ& r)  { x = mpz_get_si(r.get_mpz_t()); }
00140   inline void assign( double & r, const ZZ & z ) { r = z.get_d(); };
00141   
00142 
00143   inline void assign(QQ& x, const ZZ& r) { mpq_set_z(x.get_mpq_t(), r.get_mpz_t());}
00144   inline void assign( QQ & q, double d ) { q = d; };
00145   inline void assign( RR & r, double d ) { r = d; };
00146   inline void assign( double & r, const QQ & z ) { r = z.get_d(); };
00147   inline void assign( double & r, const RR & z ) { r = z.get_d(); };  
00148   inline void assign(RR& r, char* s)        { mpf_set_str(r.get_mpf_t(),s,10); }
00149   inline void assign(RR& x, const RR& r)    { mpf_set(x.get_mpf_t(),r.get_mpf_t()); }
00150   inline double convert( const QQ & q, type::As<double> ) { return q.get_d(); };
00151   inline double convert( const ZZ & z, type::As<double> ) { return z.get_d(); };
00152   inline double convert( const RR & r, type::As<double> ) { return r.get_d(); };
00153   //  template<class K>
00154   inline void assign( Interval<double, 3>& a, 
00155                       const QQ & b )
00156     {
00157       double  mn,mx;
00158       {
00159         numerics::rounding<double> rnd( numerics::rnd_up() );
00160         assign(mx,b);
00161       }
00162       {
00163         numerics::rounding<double> rnd( numerics::rnd_dw() );
00164         assign(mn,b);
00165       }
00166       
00167       if( mx < mn ) std::swap(mn,mx);
00168       QQ MN(mn);
00169       QQ MX(mx);
00170       new (&a) Interval<double,3>(mn,mx);
00171     };
00172 
00173 }
00174 
00175 __END_NAMESPACE_SYNAPS
00176 
00177 #endif

SYNAPS DOCUMENTATION
logo