src/etc/dec2flt_table.py

#!/usr/bin/env python3

"""
Generate powers of five using Daniel Lemire's ``Eisel-Lemire algorithm`` for use in
decimal to floating point conversions.

Specifically, computes and outputs (as Rust code) a table of 10^e for some
range of exponents e. The output is one array of 128 bit significands.
The base two exponents can be inferred using a logarithmic slope
of the decimal exponent. The approximations are normalized and rounded perfectly,
i.e., within 0.5 ULP of the true value.

Adapted from Daniel Lemire's fast_float ``table_generation.py``,
available here: <https://github.com/fastfloat/fast_float/blob/main/script/table_generation.py>.
"""
from __future__ import print_function
from math import ceil, floor, log
from collections import deque

HEADER = """
//! Pre-computed tables powers-of-5 for extended-precision representations.
//!
//! These tables enable fast scaling of the significant digits
//! of a float to the decimal exponent, with minimal rounding
//! errors, in a 128 or 192-bit representation.
//!
//! DO NOT MODIFY: Generated by `src/etc/dec2flt_table.py`
"""

STATIC_WARNING = """
// Use static to avoid long compile times: Rust compiler errors
// can have the entire table compiled multiple times, and then
// emit code multiple times, even if it's stripped out in
// the final binary.
"""

def main():
    min_exp = minimum_exponent(10)
    max_exp = maximum_exponent(10)
    bias = -minimum_exponent(5)

    print(HEADER.strip())
    print()
    print('pub const SMALLEST_POWER_OF_FIVE: i32 = {};'.format(min_exp))
    print('pub const LARGEST_POWER_OF_FIVE: i32 = {};'.format(max_exp))
    print('pub const N_POWERS_OF_FIVE: usize = ', end='')
    print('(LARGEST_POWER_OF_FIVE - SMALLEST_POWER_OF_FIVE + 1) as usize;')
    print()
    print_proper_powers(min_exp, max_exp, bias)


def minimum_exponent(base):
    return ceil(log(5e-324, base) - log(0xFFFFFFFFFFFFFFFF, base))


def maximum_exponent(base):
    return floor(log(1.7976931348623157e+308, base))


def print_proper_powers(min_exp, max_exp, bias):
    powers = deque()

    # Add negative exponents.
    # 2^(2b)/(5^−q) with b=64 + int(math.ceil(log2(5^−q)))
    powers = []
    for q in range(min_exp, 0):
        power5 = 5 ** -q
        z = 0
        while (1 << z) < power5:
            z += 1
        if q >= -27:
            b = z + 127
            c = 2 ** b // power5 + 1
            powers.append((c, q))
        else:
            b = 2 * z + 2 * 64
            c = 2 ** b // power5 + 1
            # truncate
            while c >= (1<<128):
                c //= 2
            powers.append((c, q))

    # Add positive exponents
    for q in range(0, max_exp + 1):
        power5 = 5 ** q
        # move the most significant bit in position
        while power5 < (1<<127):
            power5 *= 2
        # *truncate*
        while power5 >= (1<<128):
            power5 //= 2
        powers.append((power5, q))

    # Print the powers.
    print(STATIC_WARNING.strip())
    print('#[rustfmt::skip]')
    typ = '[(u64, u64); N_POWERS_OF_FIVE]'
    print('pub static POWER_OF_FIVE_128: {} = ['.format(typ))
    for c, exp in powers:
        hi = '0x{:x}'.format(c // (1 << 64))
        lo = '0x{:x}'.format(c % (1 << 64))
        value = '    ({}, {}), '.format(hi, lo)
        comment = '// {}^{}'.format(5, exp)
        print(value.ljust(46, ' ') + comment)
    print('];')


if __name__ == '__main__':
    main()