#!/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.

Ported from Rust's core library implementation, which itself is
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, log2
from fractions import Fraction
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/lemire_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))
    lo_mask = (1 << 64) - 1
    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()