math.fractional-power

import re
from typing import Optional

#: Every value is held as an integer count of 10^-18.
FIXED_SCALE = 10 ** 18

#: Intermediate values beyond 2^512 (about 10^136 after scaling) are refused.
_MAX_BITS = 512

_DECIMAL = re.compile(r"^[0-9]+(\.[0-9]{1,18})?$")


def _too_large(value: int) -> bool:
    return value.bit_length() > _MAX_BITS


def mul_fixed(a: int, b: int) -> int:
    """Product of two non-negative fixed-point values, floored."""
    return (a * b) // FIXED_SCALE


def pow_fixed_bounded(x: int, n: int, bound: Optional[int]) -> int:
    """x^n for a non-negative fixed-point x, by square-and-multiply from the
    lowest bit, flooring after every product.

    The order of operations is part of the contract: it is what makes
    TypeScript, Python and Rust agree to the last digit. Returns bound + 1 as
    soon as the result must exceed ``bound`` (x >= 1 only grows), so a
    bisection never builds an astronomically large number.
    """
    result = FIXED_SCALE
    base = x
    e = n
    while e > 0:
        if e % 2 == 1:
            result = mul_fixed(result, base)
            if bound is not None and result > bound and x >= FIXED_SCALE:
                return bound + 1
        e //= 2
        if e > 0:
            base = mul_fixed(base, base)
            if bound is not None and base > bound and x >= FIXED_SCALE:
                return bound + 1
            if _too_large(base):
                raise ValueError("fractional power result too large")
    if _too_large(result):
        raise ValueError("fractional power result too large")
    return result


def pow_fixed(x: int, n: int) -> int:
    """x^n in fixed point; n is a whole number of 0 or more."""
    return pow_fixed_bounded(x, n, None)


def root_fixed(x: int, q: int) -> int:
    """The q-th root of x in fixed point: the largest y with
    pow_fixed(y, q) <= x, found by bisection."""
    if q == 1:
        return x
    lo = 0
    # For x >= 1, (1 + (x - 1)/q)^q >= x (Bernoulli), so the root is below it.
    hi = FIXED_SCALE + (x - FIXED_SCALE) // q + 2 if x >= FIXED_SCALE else FIXED_SCALE + 1
    while hi - lo > 1:
        mid = (lo + hi) // 2
        if pow_fixed_bounded(mid, q, x) <= x:
            lo = mid
        else:
            hi = mid
    return lo


def fractional_power_fixed(x: int, p: int, q: int) -> int:
    """x^(p/q) in fixed point: the q-th root first, then the power, then the
    reciprocal if p < 0."""
    if isinstance(p, bool) or not isinstance(p, int) or p < -100000 or p > 100000:
        raise ValueError("exponentNumerator must be between -100000 and 100000, received %r" % (p,))
    if isinstance(q, bool) or not isinstance(q, int) or q < 1 or q > 100000:
        raise ValueError("exponentDenominator must be between 1 and 100000, received %r" % (q,))
    if x <= 0:
        raise ValueError("base must be greater than zero")
    root = root_fixed(x, q)
    if p >= 0:
        return pow_fixed(root, p)
    denominator = pow_fixed(root, -p)
    if denominator == 0:
        raise ValueError("fractional power result too large")
    result = (FIXED_SCALE * FIXED_SCALE) // denominator
    if _too_large(result):
        raise ValueError("fractional power result too large")
    return result


def parse_fixed(text: str) -> int:
    """Parse a non-negative decimal of at most 18 places into fixed point."""
    if not isinstance(text, str) or not _DECIMAL.fullmatch(text):
        raise ValueError('base must be a positive decimal with at most 18 places, received "%s"' % (text,))
    whole, _, fraction = text.partition(".")
    return int(whole) * FIXED_SCALE + int(fraction.ljust(18, "0"))


def fractional_power(base: str, exponent_numerator: int, exponent_denominator: int) -> str:
    """base^(exponent_numerator / exponent_denominator), to 12 decimal places.

    Computed in 18-place fixed point with a floor after every step, then
    rounded half-up to 12 places, so the last digit printed is right unless the
    true value sits within about 10^-15 of a rounding boundary.
    """
    x = parse_fixed(base)
    value = fractional_power_fixed(x, exponent_numerator, exponent_denominator)
    rounded = (value + 500000) // 1000000
    return "%d.%012d" % (rounded // 10 ** 12, rounded % 10 ** 12)