From 67539822c80c8057d850b542abc142c189bd3f05 Mon Sep 17 00:00:00 2001
From: Sam Moore <matches@ucc.asn.au>
Date: Fri, 3 Oct 2014 16:38:01 +0800
Subject: [PATCH] Gmprat logarithms now O(1) not O(N)

... Where the "1" is determined log(3), don't get too excited

So we can finally actually print values outside DBL_MIN->DBL_MAX
Now to go through and remove all the cast to doubles...
---
 src/debugscript.cpp          | 10 +++++++
 src/debugscript.h            |  1 +
 src/eye_of_the_rabbit.script |  1 +
 src/gmprat.h                 | 52 +++++++++++++++++-------------------
 4 files changed, 37 insertions(+), 27 deletions(-)

diff --git a/src/debugscript.cpp b/src/debugscript.cpp
index 9bb5ea6..3d13245 100644
--- a/src/debugscript.cpp
+++ b/src/debugscript.cpp
@@ -113,6 +113,12 @@ void DebugScript::ParseAction()
 	{
 		currentAction.type = AT_RecordPerformance;
 	}
+	else if (actionType == "debugfont")
+	{
+		currentAction.type = AT_DebugFont;
+		inp >> currentAction.textargs;
+	}
+
 }
 
 bool DebugScript::Execute(View *view, Screen *scr)
@@ -207,6 +213,10 @@ bool DebugScript::Execute(View *view, Screen *scr)
 	case AT_RecordPerformance:
 		PrintPerformance(view, scr);
 		break;
+	case AT_DebugFont:
+		scr->ShowDebugFont(currentAction.textargs == "1" || currentAction.textargs == "on");
+		currentAction.loops = 1;
+		break;
 	default:
 		Fatal("Unknown script command in queue.");
 	}
diff --git a/src/debugscript.h b/src/debugscript.h
index 9895ddc..d39807d 100644
--- a/src/debugscript.h
+++ b/src/debugscript.h
@@ -40,6 +40,7 @@ private:
 		AT_ClearPerformance,
 		AT_PrintPerformance,
 		AT_RecordPerformance,
+		AT_DebugFont,
 		AT_Quit
 	};
 
diff --git a/src/eye_of_the_rabbit.script b/src/eye_of_the_rabbit.script
index af3116d..b4b22c6 100644
--- a/src/eye_of_the_rabbit.script
+++ b/src/eye_of_the_rabbit.script
@@ -1,6 +1,7 @@
 # Test how well document scales back to original...
 gpu
 lazy
+debugfont off
 clearperf
 label start
 printperf
diff --git a/src/gmprat.h b/src/gmprat.h
index 6782734..3e275d2 100644
--- a/src/gmprat.h
+++ b/src/gmprat.h
@@ -8,6 +8,7 @@
 
 #include <gmp.h>
 #include <string>
+#include <climits>
 
 class Gmprat
 {
@@ -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);}
 		
-- 
2.20.1