1 from __future__ import division 2 # When true division is the default, get rid of this and add it to 3 # test_long.py instead. In the meantime, it's too obscure to try to 4 # trick just part of test_long into using future division. 5 6 import sys 7 import random 8 import math 9 import unittest 10 from test.test_support import run_unittest 11 12 # decorator for skipping tests on non-IEEE 754 platforms 13 requires_IEEE_754 = unittest.skipUnless( 14 float.__getformat__("double").startswith("IEEE"), 15 "test requires IEEE 754 doubles") 16 17 DBL_MAX = sys.float_info.max 18 DBL_MAX_EXP = sys.float_info.max_exp 19 DBL_MIN_EXP = sys.float_info.min_exp 20 DBL_MANT_DIG = sys.float_info.mant_dig 21 DBL_MIN_OVERFLOW = 2**DBL_MAX_EXP - 2**(DBL_MAX_EXP - DBL_MANT_DIG - 1) 22 23 # pure Python version of correctly-rounded true division 24 def truediv(a, b): 25 """Correctly-rounded true division for integers.""" 26 negative = a^b < 0 27 a, b = abs(a), abs(b) 28 29 # exceptions: division by zero, overflow 30 if not b: 31 raise ZeroDivisionError("division by zero") 32 if a >= DBL_MIN_OVERFLOW * b: 33 raise OverflowError("int/int too large to represent as a float") 34 35 # find integer d satisfying 2**(d - 1) <= a/b < 2**d 36 d = a.bit_length() - b.bit_length() 37 if d >= 0 and a >= 2**d * b or d < 0 and a * 2**-d >= b: 38 d += 1 39 40 # compute 2**-exp * a / b for suitable exp 41 exp = max(d, DBL_MIN_EXP) - DBL_MANT_DIG 42 a, b = a << max(-exp, 0), b << max(exp, 0) 43 q, r = divmod(a, b) 44 45 # round-half-to-even: fractional part is r/b, which is > 0.5 iff 46 # 2*r > b, and == 0.5 iff 2*r == b. 47 if 2*r > b or 2*r == b and q % 2 == 1: 48 q += 1 49 50 result = math.ldexp(float(q), exp) 51 return -result if negative else result 52 53 class TrueDivisionTests(unittest.TestCase): 54 def test(self): 55 huge = 1L << 40000 56 mhuge = -huge 57 self.assertEqual(huge / huge, 1.0) 58 self.assertEqual(mhuge / mhuge, 1.0) 59 self.assertEqual(huge / mhuge, -1.0) 60 self.assertEqual(mhuge / huge, -1.0) 61 self.assertEqual(1 / huge, 0.0) 62 self.assertEqual(1L / huge, 0.0) 63 self.assertEqual(1 / mhuge, 0.0) 64 self.assertEqual(1L / mhuge, 0.0) 65 self.assertEqual((666 * huge + (huge >> 1)) / huge, 666.5) 66 self.assertEqual((666 * mhuge + (mhuge >> 1)) / mhuge, 666.5) 67 self.assertEqual((666 * huge + (huge >> 1)) / mhuge, -666.5) 68 self.assertEqual((666 * mhuge + (mhuge >> 1)) / huge, -666.5) 69 self.assertEqual(huge / (huge << 1), 0.5) 70 self.assertEqual((1000000 * huge) / huge, 1000000) 71 72 namespace = {'huge': huge, 'mhuge': mhuge} 73 74 for overflow in ["float(huge)", "float(mhuge)", 75 "huge / 1", "huge / 2L", "huge / -1", "huge / -2L", 76 "mhuge / 100", "mhuge / 100L"]: 77 # If the "eval" does not happen in this module, 78 # true division is not enabled 79 with self.assertRaises(OverflowError): 80 eval(overflow, namespace) 81 82 for underflow in ["1 / huge", "2L / huge", "-1 / huge", "-2L / huge", 83 "100 / mhuge", "100L / mhuge"]: 84 result = eval(underflow, namespace) 85 self.assertEqual(result, 0.0, 'expected underflow to 0 ' 86 'from {!r}'.format(underflow)) 87 88 for zero in ["huge / 0", "huge / 0L", "mhuge / 0", "mhuge / 0L"]: 89 with self.assertRaises(ZeroDivisionError): 90 eval(zero, namespace) 91 92 def check_truediv(self, a, b, skip_small=True): 93 """Verify that the result of a/b is correctly rounded, by 94 comparing it with a pure Python implementation of correctly 95 rounded division. b should be nonzero.""" 96 97 a, b = long(a), long(b) 98 99 # skip check for small a and b: in this case, the current 100 # implementation converts the arguments to float directly and 101 # then applies a float division. This can give doubly-rounded 102 # results on x87-using machines (particularly 32-bit Linux). 103 if skip_small and max(abs(a), abs(b)) < 2**DBL_MANT_DIG: 104 return 105 106 try: 107 # use repr so that we can distinguish between -0.0 and 0.0 108 expected = repr(truediv(a, b)) 109 except OverflowError: 110 expected = 'overflow' 111 except ZeroDivisionError: 112 expected = 'zerodivision' 113 114 try: 115 got = repr(a / b) 116 except OverflowError: 117 got = 'overflow' 118 except ZeroDivisionError: 119 got = 'zerodivision' 120 121 self.assertEqual(expected, got, "Incorrectly rounded division {}/{}: " 122 "expected {}, got {}".format(a, b, expected, got)) 123 124 @requires_IEEE_754 125 def test_correctly_rounded_true_division(self): 126 # more stringent tests than those above, checking that the 127 # result of true division of ints is always correctly rounded. 128 # This test should probably be considered CPython-specific. 129 130 # Exercise all the code paths not involving Gb-sized ints. 131 # ... divisions involving zero 132 self.check_truediv(123, 0) 133 self.check_truediv(-456, 0) 134 self.check_truediv(0, 3) 135 self.check_truediv(0, -3) 136 self.check_truediv(0, 0) 137 # ... overflow or underflow by large margin 138 self.check_truediv(671 * 12345 * 2**DBL_MAX_EXP, 12345) 139 self.check_truediv(12345, 345678 * 2**(DBL_MANT_DIG - DBL_MIN_EXP)) 140 # ... a much larger or smaller than b 141 self.check_truediv(12345*2**100, 98765) 142 self.check_truediv(12345*2**30, 98765*7**81) 143 # ... a / b near a boundary: one of 1, 2**DBL_MANT_DIG, 2**DBL_MIN_EXP, 144 # 2**DBL_MAX_EXP, 2**(DBL_MIN_EXP-DBL_MANT_DIG) 145 bases = (0, DBL_MANT_DIG, DBL_MIN_EXP, 146 DBL_MAX_EXP, DBL_MIN_EXP - DBL_MANT_DIG) 147 for base in bases: 148 for exp in range(base - 15, base + 15): 149 self.check_truediv(75312*2**max(exp, 0), 69187*2**max(-exp, 0)) 150 self.check_truediv(69187*2**max(exp, 0), 75312*2**max(-exp, 0)) 151 152 # overflow corner case 153 for m in [1, 2, 7, 17, 12345, 7**100, 154 -1, -2, -5, -23, -67891, -41**50]: 155 for n in range(-10, 10): 156 self.check_truediv(m*DBL_MIN_OVERFLOW + n, m) 157 self.check_truediv(m*DBL_MIN_OVERFLOW + n, -m) 158 159 # check detection of inexactness in shifting stage 160 for n in range(250): 161 # (2**DBL_MANT_DIG+1)/(2**DBL_MANT_DIG) lies halfway 162 # between two representable floats, and would usually be 163 # rounded down under round-half-to-even. The tiniest of 164 # additions to the numerator should cause it to be rounded 165 # up instead. 166 self.check_truediv((2**DBL_MANT_DIG + 1)*12345*2**200 + 2**n, 167 2**DBL_MANT_DIG*12345) 168 169 # 1/2731 is one of the smallest division cases that's subject 170 # to double rounding on IEEE 754 machines working internally with 171 # 64-bit precision. On such machines, the next check would fail, 172 # were it not explicitly skipped in check_truediv. 173 self.check_truediv(1, 2731) 174 175 # a particularly bad case for the old algorithm: gives an 176 # error of close to 3.5 ulps. 177 self.check_truediv(295147931372582273023, 295147932265116303360) 178 for i in range(1000): 179 self.check_truediv(10**(i+1), 10**i) 180 self.check_truediv(10**i, 10**(i+1)) 181 182 # test round-half-to-even behaviour, normal result 183 for m in [1, 2, 4, 7, 8, 16, 17, 32, 12345, 7**100, 184 -1, -2, -5, -23, -67891, -41**50]: 185 for n in range(-10, 10): 186 self.check_truediv(2**DBL_MANT_DIG*m + n, m) 187 188 # test round-half-to-even, subnormal result 189 for n in range(-20, 20): 190 self.check_truediv(n, 2**1076) 191 192 # largeish random divisions: a/b where |a| <= |b| <= 193 # 2*|a|; |ans| is between 0.5 and 1.0, so error should 194 # always be bounded by 2**-54 with equality possible only 195 # if the least significant bit of q=ans*2**53 is zero. 196 for M in [10**10, 10**100, 10**1000]: 197 for i in range(1000): 198 a = random.randrange(1, M) 199 b = random.randrange(a, 2*a+1) 200 self.check_truediv(a, b) 201 self.check_truediv(-a, b) 202 self.check_truediv(a, -b) 203 self.check_truediv(-a, -b) 204 205 # and some (genuinely) random tests 206 for _ in range(10000): 207 a_bits = random.randrange(1000) 208 b_bits = random.randrange(1, 1000) 209 x = random.randrange(2**a_bits) 210 y = random.randrange(1, 2**b_bits) 211 self.check_truediv(x, y) 212 self.check_truediv(x, -y) 213 self.check_truediv(-x, y) 214 self.check_truediv(-x, -y) 215 216 217 def test_main(): 218 run_unittest(TrueDivisionTests) 219 220 if __name__ == "__main__": 221 test_main() 222
