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 import sys
  18 from math import ceil, log
  19 from fractions import Fraction
  20 from collections import namedtuple
  21
  22
  23 N = 64  # Size of the significand field in bits
  24 MIN_SIG = 2 ** (N - 1)
  25 MAX_SIG = (2 ** N) - 1
  26
  27 # Hand-rolled fp representation without arithmetic or any other operations.
  28 # The significand is normalized and always N bit, but the exponent is
  29 # unrestricted in range.
  30 Fp = namedtuple('Fp', 'sig exp')
  31
  32
  33 def algorithm_m(f, e):
  34     assert f > 0
  35     if e < 0:
  36         u = f
  37         v = 10 ** abs(e)
  38     else:
  39         u = f * 10 ** e
  40         v = 1
  41     k = 0
  42     x = u // v
  43     while True:
  44         if x < MIN_SIG:
  45             u <<= 1
  46             k -= 1
  47         elif x >= MAX_SIG:
  48             v <<= 1
  49             k += 1
  50         else:
  51             break
  52         x = u // v
  53     return ratio_to_float(u, v, k)
  54
  55
  56 def ratio_to_float(u, v, k):
  57     q, r = divmod(u, v)
  58     v_r = v - r
  59     z = Fp(q, k)
  60     if r < v_r:
  61         return z
  62     elif r > v_r:
  63         return next_float(z)
  64     elif q % 2 == 0:
  65         return z
  66     else:
  67         return next_float(z)
  68
  69
  70 def next_float(z):
  71     if z.sig == MAX_SIG:
  72         return Fp(MIN_SIG, z.exp + 1)
  73     else:
  74         return Fp(z.sig + 1, z.exp)
  75
  76
  77 def error(f, e, z):
  78     decimal = f * Fraction(10) ** e
  79     binary = z.sig * Fraction(2) ** z.exp
  80     abs_err = abs(decimal - binary)
  81     # The unit in the last place has value z.exp
  82     ulp_err = abs_err / Fraction(2) ** z.exp
  83     return float(ulp_err)
  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()