Clinger's correct decimal input in Python
Converts decimal strings to correctly rounded floats using big-integer arithmetic when fast paths fail.
- File
clinger_decimal_conversion.py, 81 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/clinger-decimal-conversion/. The authors are listed in prov/FZCTBB.jsonl.
"""Clinger's correct decimal input (1990): a floating-point fast path and AlgorithmM with big integers."""
import math
P = 53 # precision of binary64
K_MIN = -1074 # exponent of the least significant bit of the smallest subnormal
K_MAX = 1024 - P # exponent of the last bit of the largest finite number
FAST_EXP = 22 # 10^22 is the largest power of ten that is a binary64 number
POW10 = [10.0 ** i for i in range(FAST_EXP + 1)] # exact: 5^22 < 2^53
def algorithm_m(u: int, v: int) -> float:
"""RN(u / v) in binary64, ties to even, for integers u >= 0 and v > 0, by exact integer division.
Finds the scale k with 2^52 <= u / (v 2^k) < 2^53 (or k = K_MIN for subnormal
results), divides, and rounds by the remainder. Returns math.inf on overflow.
Cost: O(1) big-integer operations on integers of O(log u + log v) bits.
"""
if u == 0:
return 0.0
k = u.bit_length() - v.bit_length() - P # @line L5
if _scaled(u, v, k) >= 1 << P: # @line L6
k += 1 # now 2^52 <= floor(u / (v 2^k)) < 2^53
k = max(k, K_MIN) # @line L7
num, den = (u, v << k) if k >= 0 else (u << -k, v)
q, r = divmod(num, den) # @line L8
if 2 * r > den or (2 * r == den and q & 1): # @line L9
q += 1
if q == 1 << P: # rounding carried into the next binade
q, k = q >> 1, k + 1
if k > K_MAX: # @line L10
return math.inf
return math.ldexp(q, k)
def _scaled(u: int, v: int, k: int) -> int:
"""floor(u / (v 2^k)) for integers u >= 0, v > 0 and any integer k. Cost: one big-integer division."""
return (u // (v << k)) if k >= 0 else ((u << -k) // v)
def clinger_read(f: int, e: int) -> float:
"""The binary64 number nearest to f * 10^e (ties to even), for integers f >= 0 and e.
Uses the fast path when f < 2^53 and |e| <= 22, where both operands are exact
binary64 numbers and one IEEE operation rounds correctly; otherwise AlgorithmM.
Cost: O(1) for the fast path; O(1) big-integer operations on numbers of
O(log f + |e|) bits otherwise.
"""
if f < 0:
raise ValueError("f must be nonnegative; apply the sign afterwards")
if f == 0: # @line L1
return 0.0
if f < 1 << P and -FAST_EXP <= e <= FAST_EXP: # @line L2
x = float(f) # exact
return x * POW10[e] if e >= 0 else x / POW10[-e] # @line L3
if e >= 0: # @line L4
return algorithm_m(f * 10 ** e, 1)
return algorithm_m(f, 10 ** -e)
def parse_decimal(s: str) -> tuple[int, int, bool]:
"""Split a decimal string such as '-12.5e-3' into (f, e, negative) with value f * 10^e. Cost: O(len(s))."""
s = s.strip().lower()
negative = s.startswith("-")
s = s.lstrip("+-")
mantissa, _, exponent = s.partition("e")
whole, _, frac = mantissa.partition(".")
digits = (whole + frac) or "0"
if not digits.isdigit():
raise ValueError(f"not a decimal number: {s!r}")
return int(digits), (int(exponent) if exponent else 0) - len(frac), negative
def read_float(s: str) -> float:
"""Correctly rounded conversion of a decimal string to binary64. Cost: as clinger_read."""
f, e, negative = parse_decimal(s)
x = clinger_read(f, e)
return -x if negative else x
test_clinger_decimal_conversion.py 112 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/clinger-decimal-conversion/. The authors are listed in prov/FZCTBB.jsonl.
"""Tests for Clinger's decimal input: compare with an exact rational oracle and with float()."""
import math
import random
import struct
from fractions import Fraction
from clinger_decimal_conversion import algorithm_m, clinger_read, parse_decimal, read_float
MAX = struct.unpack(">d", bytes.fromhex("7fefffffffffffff"))[0]
def oracle(f: int, e: int) -> float:
"""RN(f * 10^e) through Fraction; CPython's int true division is correctly rounded."""
q = Fraction(f) * Fraction(10) ** e
try:
return q.numerator / q.denominator
except OverflowError:
return math.inf
def test_examples() -> None:
cases = {
"0.1": "0x1.999999999999ap-4",
"3.14159": "0x1.921f9f01b866ep+1",
"1e23": "0x1.52d02c7e14af6p+76",
"123456789012345678901234567890": "0x1.8ee90ff6c373ep+96",
"8.98846567431158e307": "0x1.0000000000000p+1023",
}
for s, h in cases.items():
assert read_float(s) == float.fromhex(h)
def test_ties_to_even() -> None:
assert read_float("9007199254740993") == 2.0 ** 53 # 2^53 + 1, a tie, goes down to even
assert read_float("9007199254740995") == 2.0 ** 53 + 4 # 2^53 + 3, a tie, goes up to even
assert read_float("9007199254740993.0000000000000001") == 2.0 ** 53 + 2 # just above the tie
def test_extremes() -> None:
assert read_float("1.7976931348623157e308") == MAX
assert read_float("1.7976931348623158e308") == MAX
assert read_float("1.7976931348623159e308") == math.inf
assert read_float("1e400") == math.inf
assert read_float("2.2250738585072014e-308") == 2.0 ** -1022
assert read_float("2.2250738585072011e-308") == math.nextafter(2.0 ** -1022, 0.0)
assert read_float("4.9406564584124654e-324") == 2.0 ** -1074
assert read_float("2.4703282292062327e-324") == 0.0 # below half the smallest subnormal
assert read_float("2.4703282292062328e-324") == 2.0 ** -1074 # above it
assert read_float("1e-400") == 0.0
def test_zero_and_sign() -> None:
assert read_float("0") == 0.0
assert math.copysign(1.0, read_float("-0.0")) == -1.0
assert read_float("-2.5") == -2.5
assert read_float("+1e2") == 100.0
def test_parse() -> None:
assert parse_decimal("12.50e-3") == (1250, -5, False)
assert parse_decimal("-7") == (7, 0, True)
assert parse_decimal(".5") == (5, -1, False)
def test_fast_path_matches_algorithm_m() -> None:
rng = random.Random(1)
for _ in range(3000):
f = rng.randrange(1, 1 << 53)
e = rng.randint(-22, 22)
x = clinger_read(f, e)
slow = algorithm_m(f * 10 ** e, 1) if e >= 0 else algorithm_m(f, 10 ** -e)
assert x == slow == oracle(f, e)
def test_naive_fast_path_beyond_22_is_wrong() -> None:
# Multiplying by an inexact power of ten rounds twice.
assert 3.0 * 1e23 != oracle(3, 23) # 1e23 is already rounded
assert clinger_read(3, 23) == oracle(3, 23) == 3e23
def test_random_against_oracle_and_float() -> None:
rng = random.Random(20240612)
for _ in range(4000):
digits = rng.randint(1, 40)
f = rng.randrange(10 ** (digits - 1), 10 ** digits)
e = rng.randint(-360, 300)
x = clinger_read(f, e)
assert x == oracle(f, e)
assert x == float(f"{f}e{e}")
def test_halfway_cases_from_floats() -> None:
# Exact decimal expansions of midpoints between neighbors, and values just off them.
rng = random.Random(5)
for _ in range(400):
bits = rng.getrandbits(63)
x = struct.unpack(">d", struct.pack(">Q", bits))[0]
if not math.isfinite(x) or x == 0.0 or x == MAX:
continue
y = math.nextafter(x, math.inf)
mid = (Fraction(x) + Fraction(y)) / 2
den = mid.denominator # a power of two, so mid has a finite decimal expansion
k = den.bit_length() - 1
f, e = mid.numerator * 5 ** k, -k
expected = x if struct.unpack(">Q", struct.pack(">d", x))[0] % 2 == 0 else y
assert clinger_read(f, e) == expected
assert clinger_read(f * 10 + 1, e - 1) == y
assert clinger_read(f * 10 - 1, e - 1) == x
What the next level needs
- the tests did not run in this build