[gcalctool/vala] Fix ln e

[Robert Ancell <rancell src gnome org>]

commit 2eb2c6c6280de0b3f3bfce16f96c113678d9f279
Author: Robert Ancell <robert ancell canonical com>
Date:   Sat Oct 13 22:19:03 2012 +1300

    Fix ln e

 src/number.vala |   68 +++++++++++++++++++++++++++++++++++++++++++++++-------
 1 files changed, 59 insertions(+), 9 deletions(-)
---
diff --git a/src/number.vala b/src/number.vala
index f7929fc..b8ad8aa 100644
--- a/src/number.vala
+++ b/src/number.vala
@@ -198,7 +198,7 @@ public class Number
             if (ie < 0)
             {
                 var k = -ie;
-                for (var i = 1; i <= k; ++i)
+                for (var i = 1; i <= k; i++)
                 {
                     tp <<= 4;
                     if (tp <= ib && tp != BASE && i < k)
@@ -209,7 +209,7 @@ public class Number
             }
             else if (ie > 0)
             {
-                for (var i = 1; i <= ie; ++i)
+                for (var i = 1; i <= ie; i++)
                 {
                     tp <<= 4;
                     if (tp <= ib && tp != BASE && i < ie)
@@ -2009,18 +2009,15 @@ public class Number
             /* COMPUTE FINAL CORRECTION ACCURATELY USING LNS */
             var t2 = t1.add (new Number.integer (-1));
             if (t2.is_zero () || t2.re_exponent + 1 <= 0)
-            {
-                t2 = t2.lns ();
-                return z.add (t2);
-            }
+                return z.add (t2.lns ());
 
             /* REMOVE EXPONENT TO AVOID FLOATING-POINT OVERFLOW */
             var e = t1.re_exponent;
             t1.re_exponent = 0;
-            var rx = t1.to_double ();
+            var rx = t1.to_float_old ();
             t1.re_exponent = e;
-            var rlx = Math.log (rx) + e * Math.log (BASE);
-            t2 = new Number.double (-rlx);
+            var rlx = (float) (Math.log (rx) + e * Math.log (BASE));
+            t2 = new Number.float (-(float)rlx);
 
             /* UPDATE Z AND COMPUTE ACCURATE EXP OF APPROXIMATE LOG */
             z = z.subtract (t2);
@@ -2034,6 +2031,59 @@ public class Number
         return z;
     }
 
+    // FIXME: This is here becase ln e breaks if we use the symmetric to_float
+    private float to_float_old ()
+    {
+        if (is_zero ())
+            return 0f;
+
+        var z = 0f;
+        var i = 0;
+        for (; i < T; i++)
+        {
+            z = BASE * z + re_fraction[i];
+
+            /* CHECK IF FULL SINGLE-PRECISION ACCURACY ATTAINED */
+            if (z + 1.0f <= z)
+                break;
+        }
+
+        /* NOW ALLOW FOR EXPONENT */
+        z = (float) (z * mppow_ri (BASE, re_exponent - i - 1));
+
+        if (re_sign < 0)
+            return -z;
+        else
+            return z;
+    }
+
+    private double mppow_ri (float ap, int bp)
+    {
+        if (bp == 0)
+            return 1.0f;
+
+        if (bp < 0)
+        {
+            if (ap == 0)
+                return 1.0f;
+            bp = -bp;
+            ap = 1 / ap;
+        }
+
+        var pow = 1.0;
+        while (true)
+        {
+            if ((bp & 01) != 0)
+                pow *= ap;
+            if ((bp >>= 1) != 0)
+                ap *= ap;
+            else
+                break;
+        }
+
+        return pow;
+    }
+
     /*  RETURNS MP Y = Ln (1+X) IF X IS AN MP NUMBER SATISFYING THE
      *  CONDITION ABS (X) < 1/B, ERROR OTHERWISE.
      *  USES NEWTONS METHOD TO SOLVE THE EQUATION