1 #!/usr/bin/env python2.7
4 Generate powers of ten using William Clinger's ``AlgorithmM`` for use in
5 decimal to floating point conversions.
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.
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).
16 from __future__ import print_function
18 from math import ceil, log
19 from fractions import Fraction
20 from collections import namedtuple
23 N = 64 # Size of the significand field in bits
24 MIN_SIG = 2 ** (N - 1)
25 MAX_SIG = (2 ** N) - 1
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')
33 def algorithm_m(f, e):
53 return ratio_to_float(u, v, k)
56 def ratio_to_float(u, v, k):
72 return Fp(MIN_SIG, z.exp + 1)
74 return Fp(z.sig + 1, z.exp)
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
86 //! Tables of approximations of powers of ten.
87 //! DO NOT MODIFY: Generated by `src/etc/dec2flt_table.py`
96 print_short_powers(32, 24)
98 print_short_powers(64, 53)
101 def print_proper_powers():
104 e_range = range(MIN_E, MAX_E+1)
107 z = algorithm_m(1, e)
111 print("pub const MIN_E: i16 = {};".format(MIN_E))
112 print("pub const MAX_E: i16 = {};".format(MAX_E))
114 typ = "([u64; {0}], [i16; {0}])".format(len(powers))
115 print("pub const POWERS: ", typ, " = ([", sep='')
117 print(" 0x{:x},".format(z.sig))
120 print(" {},".format(z.exp))
124 def print_short_powers(num_bits, significand_size):
125 max_sig = 2**significand_size - 1
126 # The fast path bails out for exponents >= ceil(log5(max_sig))
127 max_e = int(ceil(log(max_sig, 5)))
128 e_range = range(max_e)
129 typ = "[f{}; {}]".format(num_bits, len(e_range))
130 print("pub const F", num_bits, "_SHORT_POWERS: ", typ, " = [", sep='')
132 print(" 1e{},".format(e))
136 if __name__ == '__main__':