src/etc/dec2flt_table.py

   1 #!/usr/bin/env python2.7
   2
   3 """
   4 Generate powers of ten using William Clinger's ``AlgorithmM`` 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 64 bit significands and
   9 another array of corresponding base two exponents. The approximations are
  10 normalized and rounded perfectly, i.e., within 0.5 ULP of the true value.
  11
  12 The representation ([u64], [i16]) instead of the more natural [(u64, i16)]
  13 is used because (u64, i16) has a ton of padding which would make the table
  14 even larger, and it's already uncomfortably large (6 KiB).
  15 """
  16 from __future__ import print_function
  17 from math import ceil, log
  18 from fractions import Fraction
  19 from collections import namedtuple
  20
  21
  22 N = 64  # Size of the significand field in bits
  23 MIN_SIG = 2 ** (N - 1)
  24 MAX_SIG = (2 ** N) - 1
  25
  26 # Hand-rolled fp representation without arithmetic or any other operations.
  27 # The significand is normalized and always N bit, but the exponent is
  28 # unrestricted in range.
  29 Fp = namedtuple('Fp', 'sig exp')
  30
  31
  32 def algorithm_m(f, e):
  33     assert f > 0
  34     if e < 0:
  35         u = f
  36         v = 10 ** abs(e)
  37     else:
  38         u = f * 10 ** e
  39         v = 1
  40     k = 0
  41     x = u // v
  42     while True:
  43         if x < MIN_SIG:
  44             u <<= 1
  45             k -= 1
  46         elif x >= MAX_SIG:
  47             v <<= 1
  48             k += 1
  49         else:
  50             break
  51         x = u // v
  52     return ratio_to_float(u, v, k)
  53
  54
  55 def ratio_to_float(u, v, k):
  56     q, r = divmod(u, v)
  57     v_r = v - r
  58     z = Fp(q, k)
  59     if r < v_r:
  60         return z
  61     elif r > v_r:
  62         return next_float(z)
  63     elif q % 2 == 0:
  64         return z
  65     else:
  66         return next_float(z)
  67
  68
  69 def next_float(z):
  70     if z.sig == MAX_SIG:
  71         return Fp(MIN_SIG, z.exp + 1)
  72     else:
  73         return Fp(z.sig + 1, z.exp)
  74
  75
  76 def error(f, e, z):
  77     decimal = f * Fraction(10) ** e
  78     binary = z.sig * Fraction(2) ** z.exp
  79     abs_err = abs(decimal - binary)
  80     # The unit in the last place has value z.exp
  81     ulp_err = abs_err / Fraction(2) ** z.exp
  82     return float(ulp_err)
  83
  84
  85 HEADER = """
  86 //! Tables of approximations of powers of ten.
  87 //! DO NOT MODIFY: Generated by `src/etc/dec2flt_table.py`
  88 """
  89
  90
  91 def main():
  92     print(HEADER.strip())
  93     print()
  94     print_proper_powers()
  95     print()
  96     print_short_powers(32, 24)
  97     print()
  98     print_short_powers(64, 53)
  99
 100
 101 def print_proper_powers():
 102     MIN_E = -305
 103     MAX_E = 305
 104     e_range = range(MIN_E, MAX_E+1)
 105     powers = []
 106     for e in e_range:
 107         z = algorithm_m(1, e)
 108         err = error(1, e, z)
 109         assert err < 0.5
 110         powers.append(z)
 111     print("pub const MIN_E: i16 = {};".format(MIN_E))
 112     print("pub const MAX_E: i16 = {};".format(MAX_E))
 113     print()
 114     print("#[rustfmt::skip]")
 115     typ = "([u64; {0}], [i16; {0}])".format(len(powers))
 116     print("pub const POWERS: ", typ, " = (", sep='')
 117     print("    [")
 118     for z in powers:
 119         print("        0x{:x},".format(z.sig))
 120     print("    ],")
 121     print("    [")
 122     for z in powers:
 123         print("        {},".format(z.exp))
 124     print("    ],")
 125     print(");")
 126
 127
 128 def print_short_powers(num_bits, significand_size):
 129     max_sig = 2**significand_size - 1
 130     # The fast path bails out for exponents >= ceil(log5(max_sig))
 131     max_e = int(ceil(log(max_sig, 5)))
 132     e_range = range(max_e)
 133     typ = "[f{}; {}]".format(num_bits, len(e_range))
 134     print("#[rustfmt::skip]")
 135     print("pub const F", num_bits, "_SHORT_POWERS: ", typ, " = [", sep='')
 136     for e in e_range:
 137         print("    1e{},".format(e))
 138     print("];")
 139
 140
 141 if __name__ == '__main__':
 142     main()