[Jython-checkins] jython: Pure python version of fsum for math module.

[frank.wierzbicki]

http://hg.python.org/jython/rev/ee12ca8ebf01
changeset:   6660:ee12ca8ebf01
user:        Miki Tebeka <miki.tebeka at gmail.com>
date:        Fri May 18 14:37:25 2012 -0700
summary:
  Pure python version of fsum for math module.

files:
  Lib/_fsum.py                     |  36 ++++++++++++++++++++
  src/org/python/modules/math.java |  11 +++--
  2 files changed, 43 insertions(+), 4 deletions(-)


diff --git a/Lib/_fsum.py b/Lib/_fsum.py
new file mode 100644
--- /dev/null
+++ b/Lib/_fsum.py
@@ -0,0 +1,36 @@
+#!/usr/bin/env python
+
+from sys import float_info
+import math
+
+
+mant_dig = float_info.mant_dig
+etiny = float_info.min_exp - mant_dig
+
+def fsum(iterable):
+    """Full precision summation.  Compute sum(iterable) without any
+    intermediate accumulation of error.  Based on the 'lsum' function
+    at http://code.activestate.com/recipes/393090/
+
+    """
+    tmant, texp = 0, 0
+    for x in iterable:
+        mant, exp = math.frexp(x)
+        mant, exp = int(math.ldexp(mant, mant_dig)), exp - mant_dig
+        if texp > exp:
+            tmant <<= texp-exp
+            texp = exp
+        else:
+            mant <<= exp-texp
+        tmant += mant
+
+    # Round tmant * 2**texp to a float.  The original recipe
+    # used float(str(tmant)) * 2.0**texp for this, but that's
+    # a little unsafe because str -> float conversion can't be
+    # relied upon to do correct rounding on all platforms.
+    tail = max(len(bin(abs(tmant)))-2 - mant_dig, etiny - texp)
+    if tail > 0:
+        h = 1 << (tail-1)
+        tmant = tmant // (2*h) + bool(tmant & h and tmant & 3*h-1)
+        texp += tail
+    return math.ldexp(tmant, texp)
diff --git a/src/org/python/modules/math.java b/src/org/python/modules/math.java
--- a/src/org/python/modules/math.java
+++ b/src/org/python/modules/math.java
@@ -9,7 +9,9 @@
 import org.python.core.PyInteger;
 import org.python.core.PyLong;
 import org.python.core.PyObject;
+import org.python.core.PySystemState;
 import org.python.core.PyTuple;
+import org.python.core.__builtin__;
 import org.python.modules.math_erf;
 
 public class math implements ClassDictInit {
@@ -507,12 +509,13 @@
         return log(ONE + v);
     }
     
-    public static double fsum(double... values) {
-        // TODO need an iterable
-        // TODO need sys.float_info
-        return 0;
+    public static double fsum(final PyObject iterable) {
+        PyFloat result = (PyFloat)__builtin__.__import__("_fsum")
+            .invoke("fsum", iterable);
+        return result.asDouble();
     }
 
+
     private static double calculateLongLog(PyLong v) {
         int exp[] = new int[1];
         double x = v.scaledDoubleValue(exp);

-- 
Repository URL: http://hg.python.org/jython