]>
Commit | Line | Data |
---|---|---|
ba9703b0 | 1 | #!/usr/bin/env python3 |
e9174d1e SL |
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 | |
9cc50fc6 | 17 | from math import ceil, log |
e9174d1e SL |
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 | ||
e9174d1e SL |
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 | ||
74b04a01 | 84 | |
9cc50fc6 | 85 | HEADER = """ |
9cc50fc6 SL |
86 | //! Tables of approximations of powers of ten. |
87 | //! DO NOT MODIFY: Generated by `src/etc/dec2flt_table.py` | |
e9174d1e SL |
88 | """ |
89 | ||
9cc50fc6 | 90 | |
e9174d1e | 91 | def main(): |
9cc50fc6 SL |
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(): | |
e9174d1e SL |
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) | |
e9174d1e SL |
111 | print("pub const MIN_E: i16 = {};".format(MIN_E)) |
112 | print("pub const MAX_E: i16 = {};".format(MAX_E)) | |
113 | print() | |
60c5eb7d | 114 | print("#[rustfmt::skip]") |
9cc50fc6 | 115 | typ = "([u64; {0}], [i16; {0}])".format(len(powers)) |
6a06907d | 116 | print("pub static POWERS: ", typ, " = (", sep='') |
60c5eb7d | 117 | print(" [") |
e9174d1e | 118 | for z in powers: |
60c5eb7d XL |
119 | print(" 0x{:x},".format(z.sig)) |
120 | print(" ],") | |
121 | print(" [") | |
e9174d1e | 122 | for z in powers: |
60c5eb7d XL |
123 | print(" {},".format(z.exp)) |
124 | print(" ],") | |
125 | print(");") | |
e9174d1e SL |
126 | |
127 | ||
9cc50fc6 SL |
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)) | |
60c5eb7d | 134 | print("#[rustfmt::skip]") |
9cc50fc6 SL |
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 | ||
e9174d1e SL |
141 | if __name__ == '__main__': |
142 | main() |