• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1#!/usr/bin/env python3
2
3"""
4Generate powers of a given radix for the Bellerophon algorithm.
5
6Specifically, computes and outputs (as Rust code) a table of 10^e for some
7range of exponents e. The output is one array of 128 bit significands.
8The base two exponents can be inferred using a logarithmic slope
9of the decimal exponent. The approximations are normalized and rounded perfectly,
10i.e., within 0.5 ULP of the true value.
11
12Ported from Rust's core library implementation, which itself is
13adapted from Daniel Lemire's fast_float ``table_generation.py``,
14available here: <https://github.com/fastfloat/fast_float/blob/main/script/table_generation.py>.
15"""
16
17import math
18from collections import deque
19
20STEP_STR = "static BASE{0}_STEP: i32 = {1};"
21SMALL_MANTISSA_STR = "static BASE{0}_SMALL_MANTISSA: [u64; {1}] = ["
22SMALL_EXPONENT_STR = "static BASE{0}_SMALL_EXPONENT: [i32; {1}] = ["
23LARGE_MANTISSA_STR = "static BASE{0}_LARGE_MANTISSA: [u64; {1}] = ["
24LARGE_EXPONENT_STR = "static BASE{0}_LARGE_EXPONENT: [i32; {1}] = ["
25SMALL_INT_STR = "static BASE{0}_SMALL_INT_POWERS: [u64; {1}] = {2};"
26BIAS_STR = "static BASE{0}_BIAS: i32 = {1};"
27EXP_STR = "// {}^{}"
28POWER_STR = """pub static BASE{0}_POWERS: BellerophonPowers = BellerophonPowers {{
29    small: ExtendedFloatArray {{ mant: &BASE{0}_SMALL_MANTISSA, exp: &BASE{0}_SMALL_EXPONENT }},
30    large: ExtendedFloatArray {{ mant: &BASE{0}_LARGE_MANTISSA, exp: &BASE{0}_LARGE_EXPONENT }},
31    small_int: &BASE{0}_SMALL_INT_POWERS,
32    step: BASE{0}_STEP,
33    bias: BASE{0}_BIAS,
34}};\n"""
35
36def calculate_bitshift(base, exponent):
37    '''
38    Calculate the bitshift required for a given base. The exponent
39    is the absolute value of the max exponent (log distance from 1.)
40    '''
41
42    return 63 + math.ceil(math.log2(base**exponent))
43
44
45def next_fp(fp, base, step = 1):
46    '''Generate the next extended-floating point value.'''
47
48    return (fp[0] * (base**step), fp[1])
49
50
51def prev_fp(fp, base, step = 1):
52    '''Generate the previous extended-floating point value.'''
53
54    return (fp[0] // (base**step), fp[1])
55
56
57def normalize_fp(fp):
58    '''Normalize a extended-float so the MSB is the 64th bit'''
59
60    while fp[0] >> 64 != 0:
61        fp = (fp[0] >> 1, fp[1] + 1)
62    return fp
63
64
65def generate_small(base, count):
66    '''Generate the small powers for a given base'''
67
68    bitshift = calculate_bitshift(base, count)
69    fps = []
70    fp = (1 << bitshift, -bitshift)
71    for exp in range(count):
72        fps.append((normalize_fp(fp), exp))
73        fp = next_fp(fp, base)
74
75    # Print the small powers as integers.
76    ints = [base**i for _, i in fps]
77
78    return fps, ints
79
80
81def generate_large(base, step):
82    '''Generate the large powers for a given base.'''
83
84    # Get our starting parameters
85    min_exp = math.floor(math.log(5e-324, base) - math.log(0xFFFFFFFFFFFFFFFF, base))
86    max_exp = math.ceil(math.log(1.7976931348623157e+308, base))
87    bitshift = calculate_bitshift(base, abs(min_exp - step))
88    fps = deque()
89
90    # Add negative exponents
91    # We need to go below the minimum exponent, since we need
92    # all resulting exponents to be positive.
93    fp = (1 << bitshift, -bitshift)
94    for exp in range(-step, min_exp-step, -step):
95        fp = prev_fp(fp, base, step)
96        fps.appendleft((normalize_fp(fp), exp))
97
98    # Add positive exponents
99    fp = (1 << bitshift, -bitshift)
100    fps.append((normalize_fp(fp), 0))
101    for exp in range(step, max_exp, step):
102        fp = next_fp(fp, base, step)
103        fps.append((normalize_fp(fp), exp))
104
105    # Return the smallest exp, AKA, the bias
106    return fps, -fps[0][1]
107
108
109def print_array(base, string, fps, index):
110    '''Print an entire array'''
111
112    print(string.format(base, len(fps)))
113    for fp, exp in fps:
114        value = "    {},".format(fp[index])
115        exp = EXP_STR.format(base, exp)
116        print(value.ljust(30, " ") + exp)
117    print("];")
118
119
120def generate_base(base):
121    '''Generate all powers and variables.'''
122
123    step = math.floor(math.log(1e10, base))
124    small, ints = generate_small(base, step)
125    large, bias = generate_large(base, step)
126
127    print_array(base, SMALL_MANTISSA_STR, small, 0)
128    print_array(base, SMALL_EXPONENT_STR, small, 1)
129    print_array(base, LARGE_MANTISSA_STR, large, 0)
130    print_array(base, LARGE_EXPONENT_STR, large, 1)
131    print(SMALL_INT_STR.format(base, len(ints), ints))
132    print(STEP_STR.format(base, step))
133    print(BIAS_STR.format(base, bias))
134
135
136def generate():
137    '''Generate all bases.'''
138
139    print("// BASE{}\n".format(10))
140    generate_base(10)
141    print("")
142
143    print("// HIGH LEVEL\n// ----------\n")
144    print(POWER_STR.format(10))
145
146
147if __name__ == '__main__':
148    generate()
149