[pypy-svn] r5707 - pypy/trunk/src/pypy/appspace

mwh at codespeak.net mwh at codespeak.net
Tue Jul 27 13:39:31 CEST 2004


Author: mwh
Date: Tue Jul 27 13:39:30 2004
New Revision: 5707

Modified:
   pypy/trunk/src/pypy/appspace/_float_formatting.py
Log:
rewrite the float formatting stuff in Python as opposed to scheme dressed as
Python.


Modified: pypy/trunk/src/pypy/appspace/_float_formatting.py
==============================================================================
--- pypy/trunk/src/pypy/appspace/_float_formatting.py	(original)
+++ pypy/trunk/src/pypy/appspace/_float_formatting.py	Tue Jul 27 13:39:30 2004
@@ -8,8 +8,7 @@
 # Language Design and Implementation.
 
 # The paper contains scheme code which has been specialized for IEEE
-# doubles and converted into (still somewhat scheme-like) Python by
-# Michael Hudson.
+# doubles and converted into Python by Michael Hudson.
 
 # XXX unfortunately, we need the fixed-format output routines, the source
 # for which is not included in the paper... for now, just put up with
@@ -17,69 +16,80 @@
 
 # XXX should run this at interpreter level, really....
 
+def decode_float(f):
+    """decode_float(float) -> int, int
 
-## (define flonum->digits
-##   (lambda (v f e min-e p b B)
-##     (if (>= e 0)
-##         (if (not (= f (expt b (- p 1))))
-##             (let ([be (expt b e)])
-##               (scale (* f be 2) 2 be be 0 B v))
-##             (let* ([be (expt b e)] [be1 (* be b)])
-##                   (scale (* f be1 2) (* b 2) be1 be 0 B v)))
-##         (if (or (= e min-e) (not (= f (expt b (- p 1)))))
-##             (scale (* f 2) (* (expt b (- e)) 2) 1 1 0 B v)
-##             (scale (* f b 2) (* (expt b (- 1 e)) 2) b 1 0 B v)))))
-
-# I reject generality in the extreme: we're working with
-# ieee 754 64 bit doubles on any platform I care about.
-# This means b == 2, min-e = -1074, p = 53 above.  Also,
-# specialize for B = 10.
-
-# in:
-# v = f * 2**e
-# out:
-# [d0, d1, ..., dn], k
-# st 0.[d1][d2]...[dn] * 10**k is the "best" representation of v
+    Return (m, e), the mantissa and exponent respectively of
+    f (assuming f is an IEEE double), i.e. f == m * 2**e and
+    2**52 <= m < 2**53."""
+    m, e = math.frexp(f)
+    m = long(m*2.0**53)
+    e -= 53
+    return m, e
+
+def decompose(f):
+    """decompose(float) -> int, int, int, int
 
-def flonum2digits(v, f, e):
+    Return r, s, m_plus, m_minus for f, in the terms of
+    Burger and Dybvig's paper (see Table 1).
+
+    To spell this out: f = r/s, (r+m+)/s is halfway to the
+    next largest floating point number, (r-m-) halfway to
+    the next smallest."""
+    m, e = decode_float(f)
     if e >= 0:
-        if not f != 2**52:
+        if not m != 2**52:
             be = 2**e
-            return scale(f*be*2, 2, be, be, 0, v)
+            return m*be*2, 2, be, be
         else:
             be = 2**e
             be1 = 2*be
-            return scale(f*be1*2, 4, be1, be, 0, v)
+            return m*be1*2, 4, be1, be
     else:
-        if e == -1075 or f != 2**52:
-            return scale(f*2, 2*2**(-e), 1, 1, 0, v)
+        if e == -1075 or m != 2**52:
+            return m*2, 2**(-e+1), 1, 1
         else:
-            return scale(f*4, 2*2**(1-e), 2, 1, 0, v)
+            return m*4, 2**(-e+2), 2, 1
+
+
+def flonum2digits(f):
+    """flonum2digits(float) -> [int], int
+
+    Given a float f return [d1, ..., dn], k such that
 
+        0.[d1][d2]...[dn] * 10**k
 
-## (define generate
-##   (lambda (r s m+ m- B low-ok? high-ok?)
-##     (let ([q-r (quotient-remainder (* r B) s)]
-##           [m+ (* m+ B)]
-##           [m- (* m- B)])
-##       (let ([d (car q-r)]
-##             [r (cdr q-r)])
-##         (let ([tc1 ((if low-ok? <= <) r m-)]
-##               [tc2 ((if high-ok? >= >) (+ r m+) s)])
-##           (if (not tc1)
-##               (if (not tc2)
-##                   (cons d (generate r s m+ m- B low-ok? high-ok?))
-##                   (list (+ d 1)))
-##               (if (not tc2)
-##                   (list d)
-##                   (if (< (* r 2) s)
-##                       (list d)
-##                       (list (+ d 1))))))))))
+    is the shortest decimal representation that will
+    reproduce f when read in by a correctly rounding input
+    routine (under any strategy for breaking ties)."""
+    
+    assert f >= 0
+    if f == 0.0:
+        return ['0'], 1
+
+    # See decompose's docstring for what these mean.
+    r, s, m_plus, m_minus = decompose(f)
 
-# now the above is an example of a pointlessly recursive algorithm if
-# ever i saw one...
+    # k is the index, relative to the radix point of the
+    # most-significant non-zero digit of the infinite
+    # decimal expansion of f.  This calculation has the
+    # potential to underestimate by one (handled below).
+    k = long(math.ceil(math.log10(f) - 1e-10))
+
+    if k >= 0:
+        s *= 10 ** k
+    else:
+        scale = 10 ** -k
+        r *= scale
+        m_plus *= scale
+        m_minus *= scale
+
+    # Check that we got k right above.
+    if r + m_plus > s:
+        s *= 10
+        k += 1
 
-def generate(r, s, m_plus, m_minus):
+    # Generate the digits.
     rr = []
     while 1:
         d, r = divmod(r*10, s)
@@ -93,47 +103,8 @@
             rr.append(d)
         if tc1 or tc2:
             break
-    return rr
+        
+    assert max(rr) < 10
+    assert min(rr) >= 0
 
-
-## (define scale
-##   (lambda (r s m+ m- k B low-ok? high-ok? v)
-##     (let ([est (inexact->exact (ceiling (- (logB B v) 1e-10)))])
-##       (if (>= est 0)
-##           (fixup r (* s (exptt B est)) m+ m- est B low-ok? high-ok? )
-##           (let ([scale (exptt B (- est))])
-##             (fixup (* r scale) s (* m+ scale) (* m- scale)
-##                    est B low-ok? high-ok? ))))))
-
-def scale(r, s, m_plus, m_minus, k, v):
-    est = long(math.ceil(math.log(v, 10) - 1e-10))
-    if est >= 0:
-        return fixup(r, s * 10 ** est, m_plus, m_minus, est)
-    else:
-        scale = 10 ** -est
-        return fixup(r*scale, s, m_plus*scale, m_minus*scale, est)
-
-
-## (define fixup
-##   (lambda (r s m+ m- k B low-ok? high-ok? )
-##     (if ((if high-ok? >= >) (+ r m+) s) ; too low?
-##         (cons (+ k 1) (generate r (* s B) m+ m- B low-ok? high-ok? ))
-##         (cons k (generate r s m+ m- B low-ok? high-ok? )))))
-
-def fixup(r, s, m_plus, m_minus, k):
-    if r + m_plus > s:
-        return generate(r, s*10, m_plus, m_minus), k + 1
-    else:
-        return generate(r, s, m_plus, m_minus), k
-
-
-def float_digits(f):
-    assert f >= 0
-    if f == 0.0:
-        return ['0'], 1
-    m, e = math.frexp(f)
-    m = long(m*2.0**53)
-    e -= 53
-    ds, k = flonum2digits(f, m, e)
-    ds = map(str, ds)
-    return ds, k
+    return map(str, rr), k



More information about the Pypy-commit mailing list