1 files changed, 111 insertions, 0 deletions
diff --git a/src/etc/dec2flt_table.py b/src/etc/dec2flt_table.py
new file mode 100644
index 000000000..aa5188d96
--- /dev/null
+++ b/src/etc/dec2flt_table.py
@@ -0,0 +1,111 @@
+#!/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, log2
+from fractions import Fraction
+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))
+    lo_mask = (1 << 64) - 1
+    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()