1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
|
#!/usr/bin/env python3
"""
Generate powers of five using Daniel Lemire's ``Eisel-Lemire algorithm`` for use in
decimal to floating point conversions.
Specifically, computes and outputs (as Rust code) a table of 10^e for some
range of exponents e. The output is one array of 128 bit significands.
The base two exponents can be inferred using a logarithmic slope
of the decimal exponent. The approximations are normalized and rounded perfectly,
i.e., within 0.5 ULP of the true value.
Adapted from Daniel Lemire's fast_float ``table_generation.py``,
available here: <https://github.com/fastfloat/fast_float/blob/main/script/table_generation.py>.
"""
from __future__ import print_function
from math import ceil, floor, log
from collections import deque
HEADER = """
//! Pre-computed tables powers-of-5 for extended-precision representations.
//!
//! These tables enable fast scaling of the significant digits
//! of a float to the decimal exponent, with minimal rounding
//! errors, in a 128 or 192-bit representation.
//!
//! DO NOT MODIFY: Generated by `src/etc/dec2flt_table.py`
"""
STATIC_WARNING = """
// Use static to avoid long compile times: Rust compiler errors
// can have the entire table compiled multiple times, and then
// emit code multiple times, even if it's stripped out in
// the final binary.
"""
def main():
min_exp = minimum_exponent(10)
max_exp = maximum_exponent(10)
bias = -minimum_exponent(5)
print(HEADER.strip())
print()
print('pub const SMALLEST_POWER_OF_FIVE: i32 = {};'.format(min_exp))
print('pub const LARGEST_POWER_OF_FIVE: i32 = {};'.format(max_exp))
print('pub const N_POWERS_OF_FIVE: usize = ', end='')
print('(LARGEST_POWER_OF_FIVE - SMALLEST_POWER_OF_FIVE + 1) as usize;')
print()
print_proper_powers(min_exp, max_exp, bias)
def minimum_exponent(base):
return ceil(log(5e-324, base) - log(0xFFFFFFFFFFFFFFFF, base))
def maximum_exponent(base):
return floor(log(1.7976931348623157e+308, base))
def print_proper_powers(min_exp, max_exp, bias):
powers = deque()
# Add negative exponents.
# 2^(2b)/(5^−q) with b=64 + int(math.ceil(log2(5^−q)))
powers = []
for q in range(min_exp, 0):
power5 = 5 ** -q
z = 0
while (1 << z) < power5:
z += 1
if q >= -27:
b = z + 127
c = 2 ** b // power5 + 1
powers.append((c, q))
else:
b = 2 * z + 2 * 64
c = 2 ** b // power5 + 1
# truncate
while c >= (1<<128):
c //= 2
powers.append((c, q))
# Add positive exponents
for q in range(0, max_exp + 1):
power5 = 5 ** q
# move the most significant bit in position
while power5 < (1<<127):
power5 *= 2
# *truncate*
while power5 >= (1<<128):
power5 //= 2
powers.append((power5, q))
# Print the powers.
print(STATIC_WARNING.strip())
print('#[rustfmt::skip]')
typ = '[(u64, u64); N_POWERS_OF_FIVE]'
print('pub static POWER_OF_FIVE_128: {} = ['.format(typ))
for c, exp in powers:
hi = '0x{:x}'.format(c // (1 << 64))
lo = '0x{:x}'.format(c % (1 << 64))
value = ' ({}, {}), '.format(hi, lo)
comment = '// {}^{}'.format(5, exp)
print(value.ljust(46, ' ') + comment)
print('];')
if __name__ == '__main__':
main()
|
