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