src/etc/dec2flt_table.py

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