src/etc/dec2flt_table.py

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