Dragon4 (Steele-White) in Python
Exact big-integer digit generation for shortest float printing.
- File
dragon4.py, 94 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/dragon4/. The authors are listed in prov/17AMFE.jsonl.
"""Dragon4 free-format printing (Steele and White, 1990), in the formulation of Burger and Dybvig (1996).
Every quantity is an exact Python integer: no floating-point operation is used after the input
is decomposed, so the digits are the shortest that read back to x, and the nearest such digits.
"""
import math
import struct
P = 53 # precision of binary64
Q_MIN = -1074 # exponent of the least significant bit of a subnormal
def decompose(x: float) -> tuple[int, int]:
"""Return (m, q) with x = m * 2^q exactly, for finite x > 0, m < 2^53 and q >= -1074. Cost: O(1)."""
bits = struct.unpack(">Q", struct.pack(">d", x))[0]
biased, frac = (bits >> 52) & 0x7FF, bits & ((1 << 52) - 1)
if biased == 0:
return frac, Q_MIN
return frac | (1 << 52), biased - 1075
def dragon4(x: float) -> tuple[str, int]:
"""Shortest decimal digits of x: returns (digits, k) with x printed as 0.digits * 10^k.
Precondition: x is finite and x > 0. The digits are the fewest that lie in the rounding
interval of x (so float() reads them back to x, ties to even), and among those the nearest
to x, ties to an even last digit. Cost: O(n) big-integer operations on integers of
O(1100) binary digits for binary64, where n <= 17 is the number of digits produced.
"""
if not (math.isfinite(x) and x > 0):
raise ValueError("dragon4 needs a finite positive float")
m, q = decompose(x) # @line L1
inclusive = m % 2 == 0 # ties-to-even readers map the gap midpoints to x when m is even
if q >= 0: # @line L2
if m != 1 << (P - 1):
r, s, m_plus, m_minus = 2 * m << q, 2, 1 << q, 1 << q
else: # the gap below a power of two is half the gap above
r, s, m_plus, m_minus = 4 * m << q, 4, 2 << q, 1 << q
elif q == Q_MIN or m != 1 << (P - 1):
r, s, m_plus, m_minus = 2 * m, 2 << -q, 1, 1
else:
r, s, m_plus, m_minus = 4 * m, 4 << -q, 2, 1
# Now x = r/s, and the gaps to the midpoints with the neighbors are m_plus/s and m_minus/s.
e2 = q + m.bit_length() - 1 # 2^e2 <= x < 2^(e2 + 1)
k = (e2 * 30103) // 100000 - 1 # @line L3, k < log10(x), and k is at most 3 below the answer
if k >= 0:
s *= 10 ** k
else:
scale = 10 ** -k
r, m_plus, m_minus = r * scale, m_plus * scale, m_minus * scale
while r + m_plus > s or (inclusive and r + m_plus == s): # @line L4
s, k = s * 10, k + 1 # now every number in the rounding interval is below 10^k
digits = []
while True:
r, m_plus, m_minus = r * 10, m_plus * 10, m_minus * 10 # @line L5
d, r = divmod(r, s) # @line L6
low = r < m_minus or (inclusive and r == m_minus) # @line L7
high = r + m_plus > s or (inclusive and r + m_plus == s)
if not low and not high: # @line L8
digits.append(d)
continue
if low and high: # @line L9
if 2 * r > s or (2 * r == s and d % 2 == 1):
d += 1
elif high:
d += 1
digits.append(d)
return "".join(map(str, digits)), k # @line L10
def to_string(x: float) -> str:
"""Format x in the style of Python's repr, using dragon4 for the digits. Cost: as dragon4."""
if math.isnan(x):
return "nan"
if math.isinf(x):
return "inf" if x > 0 else "-inf"
if x == 0.0:
return "-0.0" if math.copysign(1.0, x) < 0 else "0.0"
sign = "-" if x < 0 else ""
digits, k = dragon4(abs(x))
n, point = len(digits), k # value = 0.digits * 10^k
if -4 < point <= 16:
if point <= 0:
return sign + "0." + "0" * -point + digits
if point >= n:
return sign + digits + "0" * (point - n) + ".0"
return sign + digits[:point] + "." + digits[point:]
exp = point - 1
mant = digits[0] + ("." + digits[1:] if n > 1 else "")
return f"{sign}{mant}e{'+' if exp >= 0 else '-'}{abs(exp):02d}"
test_dragon4.py 135 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/dragon4/. The authors are listed in prov/17AMFE.jsonl.
"""Tests for Dragon4: shortest, nearest, round trip, checked with exact rationals and with repr."""
import math
import random
import struct
from fractions import Fraction
import pytest
from dragon4 import decompose, dragon4, to_string
MAX = 1.7976931348623157e308
def from_bits(b: int) -> float:
return struct.unpack(">d", struct.pack(">Q", b))[0]
def interval(x: float) -> tuple[Fraction, Fraction, bool]:
"""The rounding interval of x > 0 under ties-to-even reading: (low, high, closed)."""
m, _ = decompose(x)
up = Fraction(math.nextafter(x, math.inf)) if x < MAX else Fraction(2) ** 1024
down = Fraction(math.nextafter(x, 0.0))
fx = Fraction(x)
return (fx + down) / 2, (fx + up) / 2, m % 2 == 0
def inside(y: Fraction, lo: Fraction, hi: Fraction, closed: bool) -> bool:
return lo <= y <= hi if closed else lo < y < hi
def check_shortest_and_nearest(x: float) -> None:
digits, k = dragon4(x)
n = len(digits)
assert digits[0] != "0" and digits[-1] != "0"
unit = Fraction(10) ** (k - n)
y = int(digits) * unit
lo, hi, closed = interval(x)
assert inside(y, lo, hi, closed)
fx = Fraction(x)
# No number with n - 1 significant digits is inside: such numbers below 10^k are multiples
# of 10 * unit, and it is enough to test the two that bracket x.
coarse = unit * 10
below = (fx // coarse) * coarse
assert not (below > 0 and inside(below, lo, hi, closed))
assert not inside(below + coarse, lo, hi, closed)
# Nearest among the n-digit candidates, ties to an even last digit.
for other in (y - unit, y + unit):
if other > 0 and inside(other, lo, hi, closed):
assert abs(other - fx) > abs(y - fx) or (abs(other - fx) == abs(y - fx) and int(digits) % 2 == 0)
def test_examples() -> None:
assert dragon4(0.1) == ("1", 0)
assert dragon4(0.3) == ("3", 0)
assert dragon4(0.1 + 0.2) == ("30000000000000004", 0)
assert dragon4(1e23) == ("1", 24)
assert dragon4(2.0 ** 63) == ("9223372036854776", 19)
assert dragon4(5e-324) == ("5", -323)
assert dragon4(MAX) == ("17976931348623157", 309)
assert dragon4(1.0) == ("1", 1)
def test_tie_between_two_shortest_candidates() -> None:
# 2^50 + 1/4 lies exactly halfway between ...624.2 and ...624.3; the even digit wins.
assert to_string(2.0 ** 50 + 0.25) == "1125899906842624.2"
assert to_string(2.0 ** 50 + 0.75) == "1125899906842624.8"
def test_power_of_two_boundary() -> None:
# Below a power of two the gap is half as wide; the interval is asymmetric.
for e in (-1022, -1000, -1, 0, 1, 52, 53, 100, 1023):
x = 2.0 ** e
check_shortest_and_nearest(x)
assert float(to_string(x)) == x
def test_extremes_and_specials() -> None:
for x in (5e-324, 2.0 ** -1022, math.nextafter(2.0 ** -1022, 0.0), MAX, 1e-300, 1e300):
check_shortest_and_nearest(x)
assert to_string(x) == repr(x)
assert to_string(0.0) == "0.0" and to_string(-0.0) == "-0.0"
assert to_string(math.inf) == "inf" and to_string(-math.inf) == "-inf"
assert to_string(math.nan) == "nan"
with pytest.raises(ValueError):
dragon4(0.0)
def test_random_bit_patterns_match_repr() -> None:
rng = random.Random(1990)
count = 0
while count < 20000:
x = from_bits(rng.getrandbits(64))
if math.isnan(x) or math.isinf(x):
continue
assert to_string(x) == repr(x)
count += 1
def test_small_subnormals() -> None:
# Here the rounding interval is wide (up to a factor 3 between its ends), so it can
# straddle a power of ten; compare with repr and also with every one-digit candidate.
for m in range(1, 3000):
x = m * 5e-324
assert to_string(x) == repr(x)
check_shortest_and_nearest(x)
digits, k = dragon4(x)
if len(digits) == 1:
lo, hi, closed = interval(x)
y, fx = int(digits) * Fraction(10) ** (k - 1), Fraction(x)
for z in [d * Fraction(10) ** (k - 2) for d in range(1, 10)]:
assert not inside(z, lo, hi, closed) or abs(z - fx) > abs(y - fx)
def test_random_shortest_nearest_exact() -> None:
rng = random.Random(7)
for _ in range(2000):
x = from_bits(rng.getrandbits(63))
if x == 0.0 or math.isinf(x) or math.isnan(x):
continue
check_shortest_and_nearest(x)
def test_short_decimals_print_back() -> None:
# Two decimals with at most 15 digits are farther apart than the width of a rounding
# interval of a normal binary64 number, so the shortest output is the decimal itself.
rng = random.Random(3)
for _ in range(3000):
n = rng.randint(1, 15)
f, e = rng.randrange(10 ** (n - 1), 10 ** n), rng.randint(-290, 290)
digits, k = dragon4(float(f"{f}e{e}"))
assert int(digits) * Fraction(10) ** (k - len(digits)) == f * Fraction(10) ** e
What the next level needs
- the tests did not run in this build