src/gmprat.h

   1 /**
   2  * @file gmpint.h
   3  * @brief Wraps to GMP mpq_t type using inlines
   4  */
   5
   6 #ifndef _GMPRAT_H
   7 #define _GMPRAT_H
   8
   9 #include <gmp.h>
  10 #include <string>
  11 #include <climits>
  12
  13 class Gmprat
  14 {
  15         public:
  16                 Gmprat(double d=0) {mpq_init(m_op); mpq_set_d(m_op, d);}
  17                 //Gmprat(int64_t p = 0, int64_t q = 1) {mpq_init(m_op); mpq_set_si(m_op, p, q);}
  18                 //Gmprat(const std::string & str, int base=10) {mpq_init_set_str(m_op, str.c_str(), base);}
  19                 Gmprat(const Gmprat & cpy) {mpq_init(m_op); mpq_set(m_op, cpy.m_op);}
  20                 virtual ~Gmprat() {mpq_clear(m_op);}
  21
  22                 //operator int64_t() const {return mpq_get_si(m_op);}
  23                 //operator uint64_t() const {return mpq_get_ui(m_op);}
  24                 operator double() const {return mpq_get_d(m_op);}
  25                 //operator float() const {return (float)ToDouble();}
  26                 double ToDouble() const {return mpq_get_d(m_op);}
  27                 std::string Str() const
  28                 {
  29                         //TODO: Make less hacky, if we care.
  30                         // Convert to scientific notation
  31                         if (Negative())
  32                         {
  33                                 return "-" + operator-().Str();
  34                         }
  35                         double p = Log10();
  36                         if (std::isinf(p))
  37                                 return "0";
  38
  39                         int P = (int)p;
  40                         double C = pow(10.0, p - P);
  41                         std::stringstream s;
  42                         s << C;
  43                         if (P != 0)
  44                                 s << "e"<< P;
  45                         //s << "("<<ToDouble()<<")";
  46                         return s.str();
  47                 }
  48
  49                 double LogE() const
  50                 {
  51                         // log(a/b) = log(a) - log(b)
  52                         // And if a is represented in base B as:
  53                         // a = a_N B^N + a_{N-1} B^{N-1} + ... + a_0
  54                         // => log(a) \approx log(a_N B^N)
  55                         // = log(a_N) + N log(B)
  56                         // where B is the base; ie number of values representable with one digit, ie: ULONG_MAX
  57
  58                         static double logB = log(ULONG_MAX); // compiler should optimise this anyway?
  59
  60                         // Undefined logs (should probably return NAN if -ve, not -INFINITY, but meh)
  61                         if (mpz_get_ui(mpq_numref(m_op)) == 0 || mpz_sgn(mpq_numref(m_op)) < 0)
  62                                 return -INFINITY;
  63
  64                         // Log of numerator
  65                         double lognum = log(mpq_numref(m_op)->_mp_d[abs(mpq_numref(m_op)->_mp_size) - 1]);
  66                         lognum += (abs(mpq_numref(m_op)->_mp_size)-1) * logB;
  67
  68                         // Subtract log of denominator, if it exists
  69                         // Note that denominator is not explicitly set to 1, this caused a lot of headache
  70                         if (abs(mpq_denref(m_op)->_mp_size) > 0)
  71                         {
  72                                 lognum -= log(mpq_denref(m_op)->_mp_d[abs(mpq_denref(m_op)->_mp_size)-1]);
  73                                 lognum -= (abs(mpq_denref(m_op)->_mp_size)-1) * logB;
  74                         }
  75                         return lognum;
  76                 }
  77
  78                 double Log10() const {return LogE() / log(10.0);}
  79
  80                 bool Negative() const {return (mpz_sgn(mpq_numref(m_op)) < 0);}
  81
  82                 Gmprat & operator=(const Gmprat & equ) {mpq_set(m_op, equ.m_op); return *this;}
  83                 Gmprat & operator=(const double & equ) {mpq_set_d(m_op, equ); return *this;}
  84                 Gmprat & operator+=(const Gmprat & add) {mpq_add(m_op, m_op, add.m_op); return *this;}
  85                 Gmprat & operator-=(const Gmprat & sub) {mpq_sub(m_op, m_op, sub.m_op); return *this;}
  86                 Gmprat & operator*=(const Gmprat & mul) {mpq_mul(m_op, m_op, mul.m_op); return *this;}
  87                 Gmprat & operator/=(const Gmprat & div) {mpq_div(m_op, m_op, div.m_op); return *this;}
  88
  89                 Gmprat operator+(const Gmprat & add) const {Gmprat a(*this); a += add; return a;}
  90                 Gmprat operator-(const Gmprat & sub) const {Gmprat a(*this); a -= sub; return a;}
  91                 Gmprat operator*(const Gmprat & mul) const {Gmprat a(*this); a *= mul; return a;}
  92                 Gmprat operator/(const Gmprat & div) const {Gmprat a(*this); a /= div; return a;}
  93                 //Gmprat operator%(const Gmprat & div) const {Gmprat a(*this); mpq_mod(a.m_op, a.m_op, div.m_op); return a;}
  94                 Gmprat operator-() const {return (Gmprat(0L)-*this);}
  95
  96
  97                 bool operator==(const Gmprat & cmp) const {return mpq_cmp(m_op, cmp.m_op) == 0;}
  98                 bool operator!=(const Gmprat & cmp) const {return mpq_cmp(m_op, cmp.m_op) != 0;}
  99                 bool operator<(const Gmprat & cmp) const {return mpq_cmp(m_op, cmp.m_op) < 0;}
 100                 bool operator>(const Gmprat & cmp) const {return mpq_cmp(m_op, cmp.m_op) > 0;}
 101                 bool operator<=(const Gmprat & cmp) const {return mpq_cmp(m_op, cmp.m_op) <= 0;}
 102                 bool operator>=(const Gmprat & cmp) const {return mpq_cmp(m_op, cmp.m_op) >= 0;}
 103
 104                 Gmprat Abs() const {Gmprat a(*this); mpq_abs(a.m_op, a.m_op); return a;}
 105
 106
 107                 size_t Size() const
 108                 {
 109                         return sizeof(uint64_t) * (mpq_numref(m_op)->_mp_alloc + mpq_denref(m_op)->_mp_alloc);
 110                 }
 111
 112         private:
 113                 friend std::ostream& operator<<(std::ostream& os, const Gmprat & fith);
 114                 mpq_t m_op;
 115 };
 116
 117 inline std::ostream & operator<<(std::ostream & os, const Gmprat & fith)
 118 {
 119         os << fith.Str();
 120         return os;
 121 }
 122
 123
 124 #if REALTYPE != 9
 125 inline std::string Str(const Gmprat & g) {return g.Str();}
 126 inline double Log10(const Gmprat & g) {return g.Log10();}
 127 inline size_t Size(const Gmprat & g) {return g.Size();}
 128 inline Gmprat Abs(const Gmprat & g) {return g.Abs();}
 129
 130 #endif
 131
 132 #endif //_GMPRAT_H
 133