David's final changes: more profiler features, fixes.

[ipdf/code.git] / src / gmprat.h

diff --git a/src/gmprat.h b/src/gmprat.h
index 6782734..e0c9849 100644 (file)
--- a/src/gmprat.h
+++ b/src/gmprat.h
@@ -8,6 +8,7 @@
 
 #include <gmp.h>
 #include <string>
+#include <climits>
 
 class Gmprat
 {
@@ -32,7 +33,7 @@ class Gmprat
                                return "-" + operator-().Str();
                        }
                        double p = Log10();
-                       if (isinf(p))
+                       if (std::isinf(p))
                                return "0";
                                
                        int P = (int)p;
@@ -45,39 +46,36 @@ class Gmprat
                        return s.str();
                }
                
-               double Log10() const
+               double LogE() const
                {
-                       mpz_t num; mpz_init(num); mpq_get_num(num, m_op);
-                       mpz_t den; mpz_init(den); mpq_get_den(den, m_op);
+                       // log(a/b) = log(a) - log(b)
+                       // And if a is represented in base B as:
+                       // a = a_N B^N + a_{N-1} B^{N-1} + ... + a_0
+                       // => log(a) \approx log(a_N B^N)
+                       // = log(a_N) + N log(B)
+                       // where B is the base; ie number of values representable with one digit, ie: ULONG_MAX
                        
-                       double lognum = 0;
-                       double logden = 0;
-                       while (mpz_sizeinbase(num, 10) > 10)
-                       {
-
-                               mpz_div_ui(num, num, 1e10);
-                               lognum += 10;
-                       }
-                       uint64_t n = mpz_get_ui(num);
-                       if (n == 0)
-                       {
-                               return -INFINITY;
-                       }
-                       lognum += log(n)/log(10.0);
-                       //Debug("%lu", mpz_get_ui(den));
-                       while (mpz_sizeinbase(den, 10) > 10)
+                       static double logB = log(ULONG_MAX); // compiler should optimise this anyway?
+                       
+                       // Undefined logs (should probably return NAN if -ve, not -INFINITY, but meh)
+                       if (mpz_get_ui(mpq_numref(m_op)) == 0 || mpz_sgn(mpq_numref(m_op)) < 0)
+                               return -INFINITY;                               
+       
+                       // Log of numerator
+                       double lognum = log(mpq_numref(m_op)->_mp_d[abs(mpq_numref(m_op)->_mp_size) - 1]);
+                       lognum += (abs(mpq_numref(m_op)->_mp_size)-1) * logB;
+                       
+                       // Subtract log of denominator, if it exists
+                       // Note that denominator is not explicitly set to 1, this caused a lot of headache
+                       if (abs(mpq_denref(m_op)->_mp_size) > 0)
                        {
-                               mpz_div_ui(den, den, 1e10);
-                               logden += 10;
+                               lognum -= log(mpq_denref(m_op)->_mp_d[abs(mpq_denref(m_op)->_mp_size)-1]);
+                               lognum -= (abs(mpq_denref(m_op)->_mp_size)-1) * logB;
                        }
-                       uint64_t d = mpz_get_ui(den);
-                       // if d is zero, its been rounded down we hope
-                       if (d != 0)
-                               logden += log(d)/log(10.0);
-                       
-                       return (lognum - logden);
+                       return lognum;
                }
                
+               double Log10() const {return LogE() / log(10.0);}
 
                bool Negative() const {return (mpz_sgn(mpq_numref(m_op)) < 0);}
                
@@ -106,6 +104,11 @@ class Gmprat
                Gmprat Abs() const {Gmprat a(*this); mpq_abs(a.m_op, a.m_op); return a;}
                
                
+               size_t Size() const
+               {
+                       return sizeof(uint64_t) * (mpq_numref(m_op)->_mp_alloc + mpq_denref(m_op)->_mp_alloc);
+               }
+               
        private:
                friend std::ostream& operator<<(std::ostream& os, const Gmprat & fith);
                mpq_t m_op;
@@ -118,7 +121,13 @@ inline std::ostream & operator<<(std::ostream & os, const Gmprat & fith)
 }
 
 
+#if REALTYPE != 9
+inline std::string Str(const Gmprat & g) {return g.Str();}
+inline double Log10(const Gmprat & g) {return g.Log10();}
+inline size_t Size(const Gmprat & g) {return g.Size();}
+inline Gmprat Abs(const Gmprat & g) {return g.Abs();}
 
+#endif
 
 #endif //_GMPRAT_H