[pypy-svn] r75746 - in pypy/branch/fast-forward/pypy: objspace/std rlib/rstruct rlib/rstruct/test rlib/test
benjamin at codespeak.net
benjamin at codespeak.net
Thu Jul 1 23:59:29 CEST 2010
Author: benjamin
Date: Thu Jul 1 23:59:27 2010
New Revision: 75746
Added:
pypy/branch/fast-forward/pypy/rlib/rstruct/test/
pypy/branch/fast-forward/pypy/rlib/rstruct/test/test_ieee.py (contents, props changed)
Removed:
pypy/branch/fast-forward/pypy/rlib/test/test_ieee.py
Modified:
pypy/branch/fast-forward/pypy/objspace/std/marshal_impl.py
pypy/branch/fast-forward/pypy/rlib/rstruct/ieee.py
Log:
correct code for packing and unpacking floats
This is courtesy Mark Dickinson, to whom we now owe much.
Modified: pypy/branch/fast-forward/pypy/objspace/std/marshal_impl.py
==============================================================================
--- pypy/branch/fast-forward/pypy/objspace/std/marshal_impl.py (original)
+++ pypy/branch/fast-forward/pypy/objspace/std/marshal_impl.py Thu Jul 1 23:59:27 2010
@@ -169,11 +169,11 @@
def pack_float(f):
result = []
- ieee.pack_float(result, f, 8, False)
+ ieee.pack_float8(result, f)
return ''.join(result)
def unpack_float(s):
- return ieee.unpack_float(s, False)
+ return ieee.unpack_float8(s)
def marshal_w__Float(space, w_float, m):
if m.version > 1:
Modified: pypy/branch/fast-forward/pypy/rlib/rstruct/ieee.py
==============================================================================
--- pypy/branch/fast-forward/pypy/rlib/rstruct/ieee.py (original)
+++ pypy/branch/fast-forward/pypy/rlib/rstruct/ieee.py Thu Jul 1 23:59:27 2010
@@ -3,107 +3,135 @@
"""
import math
-from pypy.rlib.rarithmetic import r_longlong, isinf, isnan, INFINITY, NAN
-def pack_float(result, number, size, bigendian):
- """Append to 'result' the 'size' characters of the 32-bit or 64-bit
- IEEE representation of the number.
+from pypy.rlib import rarithmetic
+
+
+def round_to_nearest(x):
+ """Python 3 style round: round a float x to the nearest int, but
+ unlike the builtin Python 2.x round function:
+
+ - return an int, not a float
+ - do round-half-to-even, not round-half-away-from-zero.
+
+ We assume that x is finite and nonnegative; except wrong results
+ if you use this for negative x.
+
"""
- if size == 4:
- bias = 127
- exp = 8
- prec = 23
- else:
- bias = 1023
- exp = 11
- prec = 52
-
- if isnan(number):
- sign = 0x80
- man, e = 1.5, bias + 1
+ int_part = int(x)
+ frac_part = x - int_part
+ if frac_part > 0.5 or frac_part == 0.5 and int_part & 1 == 1:
+ int_part += 1
+ return int_part
+
+
+def float_unpack(Q, size):
+ """Convert a 32-bit or 64-bit integer created
+ by float_pack into a Python float."""
+
+ if size == 8:
+ MIN_EXP = -1021 # = sys.float_info.min_exp
+ MAX_EXP = 1024 # = sys.float_info.max_exp
+ MANT_DIG = 53 # = sys.float_info.mant_dig
+ BITS = 64
+ elif size == 4:
+ MIN_EXP = -125 # C's FLT_MIN_EXP
+ MAX_EXP = 128 # FLT_MAX_EXP
+ MANT_DIG = 24 # FLT_MANT_DIG
+ BITS = 32
else:
- if number < 0:
- sign = 0x80
- number *= -1
- elif number == 0.0:
- for i in range(size):
- result.append('\x00')
- return
- else:
- sign = 0x00
- if isinf(number):
- man, e = 1.0, bias + 1
- else:
- man, e = math.frexp(number)
+ raise ValueError("invalid size value")
- if 0.5 <= man and man < 1.0:
- man *= 2
- e -= 1
- man -= 1
- e += bias
- power_of_two = r_longlong(1) << prec
- mantissa = r_longlong(power_of_two * man + 0.5)
- if mantissa >> prec :
- mantissa = 0
- e += 1
-
- for i in range(size-2):
- result.append(chr(mantissa & 0xff))
- mantissa >>= 8
- x = (mantissa & ((1<<(15-exp))-1)) | ((e & ((1<<(exp-7))-1))<<(15-exp))
- result.append(chr(x))
- x = sign | e >> (exp - 7)
- result.append(chr(x))
- if bigendian:
- first = len(result) - size
- last = len(result) - 1
- for i in range(size // 2):
- (result[first + i], result[last - i]) = (
- result[last - i], result[first + i])
-
-def unpack_float(input, bigendian):
- """Interpret the 'input' string into a 32-bit or 64-bit
- IEEE representation a the number.
- """
- size = len(input)
- bytes = []
- if bigendian:
- reverse_mask = size - 1
+ if Q >> BITS:
+ raise ValueError("input out of range")
+
+ # extract pieces
+ sign = Q >> BITS - 1
+ exp = (Q & ((1 << BITS - 1) - (1 << MANT_DIG - 1))) >> MANT_DIG - 1
+ mant = Q & ((1 << MANT_DIG - 1) - 1)
+
+ if exp == MAX_EXP - MIN_EXP + 2:
+ # nan or infinity
+ result = float('nan') if mant else float('inf')
+ elif exp == 0:
+ # subnormal or zero
+ result = math.ldexp(float(mant), MIN_EXP - MANT_DIG)
else:
- reverse_mask = 0
- nonzero = False
- for i in range(size):
- x = ord(input[i ^ reverse_mask])
- bytes.append(x)
- nonzero |= x
- if not nonzero:
- return 0.0
- if size == 4:
- bias = 127
- exp = 8
- prec = 23
+ # normal
+ mant += 1 << MANT_DIG - 1
+ result = math.ldexp(float(mant), exp + MIN_EXP - MANT_DIG - 1)
+ return -result if sign else result
+
+
+def float_pack(x, size):
+ """Convert a Python float x into a 64-bit unsigned integer
+ with the same byte representation."""
+
+ if size == 8:
+ MIN_EXP = -1021 # = sys.float_info.min_exp
+ MAX_EXP = 1024 # = sys.float_info.max_exp
+ MANT_DIG = 53 # = sys.float_info.mant_dig
+ BITS = 64
+ elif size == 4:
+ MIN_EXP = -125 # C's FLT_MIN_EXP
+ MAX_EXP = 128 # FLT_MAX_EXP
+ MANT_DIG = 24 # FLT_MANT_DIG
+ BITS = 32
else:
- bias = 1023
- exp = 11
- prec = 52
- mantissa_scale_factor = 0.5 ** prec # this is constant-folded if it's
- # right after the 'if'
- mantissa = r_longlong(bytes[size-2] & ((1<<(15-exp))-1))
- for i in range(size-3, -1, -1):
- mantissa = mantissa << 8 | bytes[i]
- mantissa = 1 + mantissa * mantissa_scale_factor
- mantissa *= 0.5
- e = (bytes[-1] & 0x7f) << (exp - 7)
- e += (bytes[size-2] >> (15 - exp)) & ((1<<(exp - 7)) -1)
- e -= bias
- e += 1
- sign = bytes[-1] & 0x80
- if e == bias + 2:
- if mantissa == 0.5:
- number = INFINITY
- else:
- return NAN
+ raise ValueError("invalid size value")
+
+ sign = rarithmetic.copysign(1.0, x) < 0.0
+ if rarithmetic.isinf(x):
+ mant = 0
+ exp = MAX_EXP - MIN_EXP + 2
+ elif rarithmetic.isnan(x):
+ mant = 1 << (MANT_DIG-2) # other values possible
+ exp = MAX_EXP - MIN_EXP + 2
+ elif x == 0.0:
+ mant = 0
+ exp = 0
else:
- number = math.ldexp(mantissa,e)
- if sign : number = -number
- return number
+ m, e = math.frexp(abs(x)) # abs(x) == m * 2**e
+ exp = e - (MIN_EXP - 1)
+ if exp > 0:
+ # Normal case.
+ mant = round_to_nearest(m * (1 << MANT_DIG))
+ mant -= 1 << MANT_DIG - 1
+ else:
+ # Subnormal case.
+ if exp + MANT_DIG - 1 >= 0:
+ mant = round_to_nearest(m * (1 << exp + MANT_DIG - 1))
+ else:
+ mant = 0
+ exp = 0
+
+ # Special case: rounding produced a MANT_DIG-bit mantissa.
+ assert 0 <= mant <= 1 << MANT_DIG - 1
+ if mant == 1 << MANT_DIG - 1:
+ mant = 0
+ exp += 1
+
+ # Raise on overflow (in some circumstances, may want to return
+ # infinity instead).
+ if exp >= MAX_EXP - MIN_EXP + 2:
+ raise OverflowError("float too large to pack in this format")
+
+ # check constraints
+ assert 0 <= mant < 1 << MANT_DIG - 1
+ assert 0 <= exp <= MAX_EXP - MIN_EXP + 2
+ assert 0 <= sign <= 1
+ return ((sign << BITS - 1) | (exp << MANT_DIG - 1)) | mant
+
+
+def pack_float8(result, x):
+ unsigned = float_pack(x, 8)
+ print unsigned
+ for i in range(8):
+ result.append(chr((unsigned >> (i * 8)) & 0xFF))
+
+
+def unpack_float8(s):
+ unsigned = 0
+ for i in range(8):
+ unsigned |= ord(s[i]) << (i * 8)
+ return float_unpack(unsigned, 8)
Added: pypy/branch/fast-forward/pypy/rlib/rstruct/test/test_ieee.py
==============================================================================
--- (empty file)
+++ pypy/branch/fast-forward/pypy/rlib/rstruct/test/test_ieee.py Thu Jul 1 23:59:27 2010
@@ -0,0 +1,103 @@
+import random
+import struct
+
+from pypy.rlib.rarithmetic import isnan
+from pypy.rlib.rstruct.ieee import float_pack, float_unpack
+
+
+class TestFloatPacking:
+
+ def check_float(self, x):
+ # check roundtrip
+ Q = float_pack(x, 8)
+ y = float_unpack(Q, 8)
+ assert repr(x) == repr(y)
+
+ # check that packing agrees with the struct module
+ struct_pack8 = struct.unpack('<Q', struct.pack('<d', x))[0]
+ float_pack8 = float_pack(x, 8)
+ assert struct_pack8 == float_pack8
+
+ # check that packing agrees with the struct module
+ try:
+ struct_pack4 = struct.unpack('<L', struct.pack('<f', x))[0]
+ except OverflowError:
+ struct_pack4 = "overflow"
+ try:
+ float_pack4 = float_pack(x, 4)
+ except OverflowError:
+ float_pack4 = "overflow"
+ assert struct_pack4 == float_pack4
+
+ # if we didn't overflow, try round-tripping the binary32 value
+ if float_pack4 != "overflow":
+ roundtrip = float_pack(float_unpack(float_pack4, 4), 4)
+ assert float_pack4 == roundtrip
+
+ def test_infinities(self):
+ self.check_float(float('inf'))
+ self.check_float(float('-inf'))
+
+ def test_zeros(self):
+ self.check_float(0.0)
+ self.check_float(-0.0)
+
+ def test_nans(self):
+ self.check_float(float('nan'))
+ self.check_float(float('-nan'))
+
+ def test_simple(self):
+ test_values = [1e-10, 0.00123, 0.5, 0.7, 1.0, 123.456, 1e10]
+ for value in test_values:
+ self.check_float(value)
+ self.check_float(-value)
+
+ def test_subnormal(self):
+ # special boundaries
+ self.check_float(2**-1074)
+ self.check_float(2**-1022)
+ self.check_float(2**-1021)
+ self.check_float((2**53-1)*2**-1074)
+ self.check_float((2**52-1)*2**-1074)
+ self.check_float((2**52+1)*2**-1074)
+
+ # other subnormals
+ self.check_float(1e-309)
+ self.check_float(1e-320)
+
+ def test_powers_of_two(self):
+ # exact powers of 2
+ for k in range(-1074, 1024):
+ self.check_float(2.**k)
+
+ # and values near powers of 2
+ for k in range(-1074, 1024):
+ self.check_float((2 - 2**-52) * 2.**k)
+
+ def test_float4_boundaries(self):
+ # Exercise IEEE 754 binary32 boundary cases.
+ self.check_float(2**128.)
+ # largest representable finite binary32 value
+ self.check_float((1 - 2.**-24) * 2**128.)
+ # halfway case: rounds up to an overflowing value
+ self.check_float((1 - 2.**-25) * 2**128.)
+ self.check_float(2**-125)
+ # smallest normal
+ self.check_float(2**-126)
+ # smallest positive binary32 value (subnormal)
+ self.check_float(2**-149)
+ # 2**-150 should round down to 0
+ self.check_float(2**-150)
+ # but anything even a tiny bit larger should round up to 2**-149
+ self.check_float((1 + 2**-52) * 2**-150)
+
+ def test_random(self):
+ # construct a Python float from random integer, using struct
+ for _ in xrange(100000):
+ Q = random.randrange(2**64)
+ x = struct.unpack('<d', struct.pack('<Q', Q))[0]
+ # nans are tricky: we can't hope to reproduce the bit
+ # pattern exactly, so check_float will fail for a random nan.
+ if isnan(x):
+ continue
+ self.check_float(x)
More information about the Pypy-commit
mailing list