X-Git-Url: https://git.ucc.asn.au/?p=ipdf%2Fcode.git;a=blobdiff_plain;f=src%2Fgmprat.h;h=3e275d28a61e739b7650eb3a91e4b2f7d1b20667;hp=5df43ef85f77e33e9eb3527dd959b433e665601a;hb=67539822c80c8057d850b542abc142c189bd3f05;hpb=ef0af5fd0129161a9e079bd3cd1298b53f1fe11a

diff --git a/src/gmprat.h b/src/gmprat.h
index 5df43ef..3e275d2 100644
--- a/src/gmprat.h
+++ b/src/gmprat.h
@@ -8,6 +8,7 @@
 
 #include <gmp.h>
 #include <string>
+#include <climits>
 
 class Gmprat
 {
@@ -20,17 +21,64 @@ class Gmprat
 		
 		//operator int64_t() const {return mpq_get_si(m_op);}
 		//operator uint64_t() const {return mpq_get_ui(m_op);}
-		//operator double() const {return mpq_get_d(m_op);}
+		operator double() const {return mpq_get_d(m_op);}
+		//operator float() const {return (float)ToDouble();}
 		double ToDouble() const {return mpq_get_d(m_op);}
-		std::string Str(int base = 10) const
+		std::string Str() const
 		{
 			//TODO: Make less hacky, if we care.
-			char * buff = mpq_get_str(NULL, 10, m_op);
-			std::string result(buff);
-			free(buff);
-			return result;
+			// Convert to scientific notation
+			if (Negative())
+			{
+				return "-" + operator-().Str();
+			}
+			double p = Log10();
+			if (isinf(p))
+				return "0";
+				
+			int P = (int)p;
+			double C = pow(10.0, p - P);
+			std::stringstream s;
+			s << C;
+			if (P != 0)
+				s << "e"<< P;
+			//s << "("<<ToDouble()<<")";
+			return s.str();
 		}
 		
+		double LogE() const
+		{
+			// 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
+			
+			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)
+			{
+				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;
+			}
+			return lognum;
+		}
+		
+		double Log10() const {return LogE() / log(10.0);}
+
+		bool Negative() const {return (mpz_sgn(mpq_numref(m_op)) < 0);}
+		
 		Gmprat & operator=(const Gmprat & equ) {mpq_set(m_op, equ.m_op); return *this;}
 		Gmprat & operator=(const double & equ) {mpq_set_d(m_op, equ); return *this;}
 		Gmprat & operator+=(const Gmprat & add) {mpq_add(m_op, m_op, add.m_op); return *this;}
@@ -57,9 +105,16 @@ class Gmprat
 		
 		
 	private:
+		friend std::ostream& operator<<(std::ostream& os, const Gmprat & fith);
 		mpq_t m_op;
 };	
 
+inline std::ostream & operator<<(std::ostream & os, const Gmprat & fith)
+{
+	os << fith.Str();
+	return os;
+}
+