Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# This file is part of Hypothesis, which may be found at
2# https://github.com/HypothesisWorks/hypothesis/
3#
4# Copyright the Hypothesis Authors.
5# Individual contributors are listed in AUTHORS.rst and the git log.
6#
7# This Source Code Form is subject to the terms of the Mozilla Public License,
8# v. 2.0. If a copy of the MPL was not distributed with this file, You can
9# obtain one at https://mozilla.org/MPL/2.0/.
11from array import array
13from hypothesis.internal.floats import float_to_int, int_to_float
15"""
16This module implements support for arbitrary floating point numbers in
17Conjecture. It doesn't make any attempt to get a good distribution, only to
18get a format that will shrink well.
20It works by defining an encoding of non-negative floating point numbers
21(including NaN values with a zero sign bit) that has good lexical shrinking
22properties.
24This encoding is a tagged union of two separate encodings for floating point
25numbers, with the tag being the first bit of 64 and the remaining 63-bits being
26the payload.
28If the tag bit is 0, the next 7 bits are ignored, and the remaining 7 bytes are
29interpreted as a 7 byte integer in big-endian order and then converted to a
30float (there is some redundancy here, as 7 * 8 = 56, which is larger than the
31largest integer that floating point numbers can represent exactly, so multiple
32encodings may map to the same float).
34If the tag bit is 1, we instead use something that is closer to the normal
35representation of floats (and can represent every non-negative float exactly)
36but has a better ordering:
381. NaNs are ordered after everything else.
392. Infinity is ordered after every finite number.
403. The sign is ignored unless two floating point numbers are identical in
41 absolute magnitude. In that case, the positive is ordered before the
42 negative.
434. Positive floating point numbers are ordered first by int(x) where
44 encoding(x) < encoding(y) if int(x) < int(y).
455. If int(x) == int(y) then x and y are sorted towards lower denominators of
46 their fractional parts.
48The format of this encoding of floating point goes as follows:
50 [exponent] [mantissa]
52Each of these is the same size their equivalent in IEEE floating point, but are
53in a different format.
55We translate exponents as follows:
57 1. The maximum exponent (2 ** 11 - 1) is left unchanged.
58 2. We reorder the remaining exponents so that all of the positive exponents
59 are first, in increasing order, followed by all of the negative
60 exponents in decreasing order (where positive/negative is done by the
61 unbiased exponent e - 1023).
63We translate the mantissa as follows:
65 1. If the unbiased exponent is <= 0 we reverse it bitwise.
66 2. If the unbiased exponent is >= 52 we leave it alone.
67 3. If the unbiased exponent is in the range [1, 51] then we reverse the
68 low k bits, where k is 52 - unbiased exponent.
70The low bits correspond to the fractional part of the floating point number.
71Reversing it bitwise means that we try to minimize the low bits, which kills
72off the higher powers of 2 in the fraction first.
73"""
76MAX_EXPONENT = 0x7FF
78BIAS = 1023
79MAX_POSITIVE_EXPONENT = MAX_EXPONENT - 1 - BIAS
82def exponent_key(e: int) -> float:
83 if e == MAX_EXPONENT:
84 return float("inf")
85 unbiased = e - BIAS
86 if unbiased < 0:
87 return 10000 - unbiased
88 else:
89 return unbiased
92ENCODING_TABLE = array("H", sorted(range(MAX_EXPONENT + 1), key=exponent_key))
93DECODING_TABLE = array("H", [0]) * len(ENCODING_TABLE)
95for i, b in enumerate(ENCODING_TABLE):
96 DECODING_TABLE[b] = i
98del i, b
101def decode_exponent(e: int) -> int:
102 """Take an integer and turn it into a suitable floating point exponent
103 such that lexicographically simpler leads to simpler floats."""
104 assert 0 <= e <= MAX_EXPONENT
105 return ENCODING_TABLE[e]
108def encode_exponent(e: int) -> int:
109 """Take a floating point exponent and turn it back into the equivalent
110 result from conjecture."""
111 assert 0 <= e <= MAX_EXPONENT
112 return DECODING_TABLE[e]
115def reverse_byte(b: int) -> int:
116 result = 0
117 for _ in range(8):
118 result <<= 1
119 result |= b & 1
120 b >>= 1
121 return result
124# Table mapping individual bytes to the equivalent byte with the bits of the
125# byte reversed. e.g. 1=0b1 is mapped to 0xb10000000=0x80=128. We use this
126# precalculated table to simplify calculating the bitwise reversal of a longer
127# integer.
128REVERSE_BITS_TABLE = bytearray(map(reverse_byte, range(256)))
131def reverse64(v: int) -> int:
132 """Reverse a 64-bit integer bitwise.
134 We do this by breaking it up into 8 bytes. The 64-bit integer is then the
135 concatenation of each of these bytes. We reverse it by reversing each byte
136 on its own using the REVERSE_BITS_TABLE above, and then concatenating the
137 reversed bytes.
139 In this case concatenating consists of shifting them into the right
140 position for the word and then oring the bits together.
141 """
142 assert v.bit_length() <= 64
143 return (
144 (REVERSE_BITS_TABLE[(v >> 0) & 0xFF] << 56)
145 | (REVERSE_BITS_TABLE[(v >> 8) & 0xFF] << 48)
146 | (REVERSE_BITS_TABLE[(v >> 16) & 0xFF] << 40)
147 | (REVERSE_BITS_TABLE[(v >> 24) & 0xFF] << 32)
148 | (REVERSE_BITS_TABLE[(v >> 32) & 0xFF] << 24)
149 | (REVERSE_BITS_TABLE[(v >> 40) & 0xFF] << 16)
150 | (REVERSE_BITS_TABLE[(v >> 48) & 0xFF] << 8)
151 | (REVERSE_BITS_TABLE[(v >> 56) & 0xFF] << 0)
152 )
155MANTISSA_MASK = (1 << 52) - 1
158def reverse_bits(x: int, n: int) -> int:
159 assert x.bit_length() <= n <= 64
160 x = reverse64(x)
161 x >>= 64 - n
162 return x
165def update_mantissa(unbiased_exponent: int, mantissa: int) -> int:
166 if unbiased_exponent <= 0:
167 mantissa = reverse_bits(mantissa, 52)
168 elif unbiased_exponent <= 51:
169 n_fractional_bits = 52 - unbiased_exponent
170 fractional_part = mantissa & ((1 << n_fractional_bits) - 1)
171 mantissa ^= fractional_part
172 mantissa |= reverse_bits(fractional_part, n_fractional_bits)
173 return mantissa
176def lex_to_float(i: int) -> float:
177 assert i.bit_length() <= 64
178 has_fractional_part = i >> 63
179 if has_fractional_part:
180 exponent = (i >> 52) & ((1 << 11) - 1)
181 exponent = decode_exponent(exponent)
182 mantissa = i & MANTISSA_MASK
183 mantissa = update_mantissa(exponent - BIAS, mantissa)
185 assert mantissa.bit_length() <= 52
187 return int_to_float((exponent << 52) | mantissa)
188 else:
189 integral_part = i & ((1 << 56) - 1)
190 return float(integral_part)
193def float_to_lex(f: float) -> int:
194 if is_simple(f):
195 assert f >= 0
196 return int(f)
197 return base_float_to_lex(f)
200def base_float_to_lex(f: float) -> int:
201 i = float_to_int(f)
202 i &= (1 << 63) - 1
203 exponent = i >> 52
204 mantissa = i & MANTISSA_MASK
205 mantissa = update_mantissa(exponent - BIAS, mantissa)
206 exponent = encode_exponent(exponent)
208 assert mantissa.bit_length() <= 52
209 return (1 << 63) | (exponent << 52) | mantissa
212def is_simple(f: float) -> int:
213 try:
214 i = int(f)
215 except (ValueError, OverflowError):
216 return False
217 if i != f:
218 return False
219 return i.bit_length() <= 56