I’m now thinking that math.fdiv
would make some sense, as a counterpart to math.fmod
. The math.fmod
docstring documents it as returning something “equal to x - n*y
for some integer n
”, and it would be nice to be able to actually get hold of that integer n
.
Below is some non-Fraction
-based Python code that could (after translation to C, of course) form the basis of such a function.
First, here’s an inner function that does the job of computing both quotient and remainder when dividing a float x
by another float y
under some restrictions - namely, that both x
and y
should be finite, that y
should be positive, and that x / y
shouldn’t be too close to overflowing. It guarantees that the returned remainder is smaller than y
in absolute value, but doesn’t make guarantees about the sign of that remainder:
def _fdiv_inner(x: float, y: float) -> tuple[int, float]:
"""
Divide one float by another, giving integer quotient and float remainder.
Assumes y is finite and positive, x is finite, and x / y does not overflow.
Returns a pair (q, r) such that x == q * y + r (exactly) and |r| < |y|.
"""
acc = 0
while y <= abs(x):
q = round(x / y)
# Could use a single fma in place of the two lines below, if supported:
# x := -fma(q, y, -x)
xlo, xhigh = exact_mult(q, y)
x = x - xhigh - xlo
acc += q
return acc, x
Here exact_mult
multiplies q
by y
and gives an exact answer in the form of the sum of two float
s: a low and a high part; the code for that is below. We could use a single fused-multiply-add for this, where available. The key point here (which requires some analysis and proof), is that after q = round(x / y)
, the value x - q * y
is always exactly representable as a float
with no precision loss; then all we need to do is compute it.
Given this restricted inner function, here’s an outer function that drives it, deals with special cases, and the like. In particular, we have to do some work to deal with the possibility that x / y
might overflow.
def fdiv(x: float, y: float) -> tuple[int, float]:
"""
Divide one float by another, giving integer quotient and float remainder.
The result (q, r) satisfies x == q * y + r (exactly), and abs(r) < abs(y),
with r having the same sign as x.
Raises ValueError if either x or y is not finite, and ZeroDivisionError
if x is finite and y is zero.
"""
# Handle error cases.
if not (math.isfinite(x) and math.isfinite(y)):
raise ValueError("x and y must be finite")
if not y:
raise ZeroDivisionError("Division by zero")
# Reduce to the case x and y finite, y positive, x nonnegative.
sign_x, x = sign(x), abs(x)
sign_y, y = sign(y), abs(y)
# If x is *much* larger than y, x / y will overflow. In that case we do a
# form of long division: first divide with a shifted version of x, then
# shift back and continue with the remainder. In the common case where the
# exponents of x and y are not too far apart, the for loop body is executed
# only once. In extreme cases, it'll be executed two or three times.
exp_diff = math.frexp(x)[1] - math.frexp(y)[1]
x_shifts = list(range(exp_diff - 954, 0, -954)) + [0]
acc = 0
for x_shift in x_shifts:
q, r = _fdiv_inner(math.ldexp(x, -x_shift), y)
x = math.ldexp(r, x_shift)
acc += q << x_shift
# Ensure that the returned remainder has the same sign as the original x.
if x < 0.0:
x += y
acc -= 1
x += 0.0
return acc * sign_x * sign_y, x * sign_x
Here sign
does the obvious thing, returning the integer -1
if its argument has its sign bit set (so including the case of negative zero) and 1 otherwise.
Here’s the whole thing, including the various helper functions.
import math
C = float.fromhex("0x1.0000002000000p+27")
def sign(x: float) -> int:
return int(math.copysign(1.0, x))
def split(x: float) -> tuple[float, float]:
"""Veltkamp splitting. Returns low, high parts."""
p = C * x
h = p + (x - p)
return x - h, h
def exact_mult(x: float, y: float) -> tuple[float, float]:
"""Dekker exact multiplication. Returns low, high parts of the product."""
xl, xh = split(x)
yl, yh = split(y)
h = x * y
return -h + xh * yh + xh * yl + xl * yh + xl * yl, h
def fdiv(x: float, y: float) -> tuple[int, float]:
"""
Divide one float by another, giving integer quotient and float remainder.
The result (q, r) satisfies x == q * y + r (exactly), and abs(r) < abs(y),
with r having the same sign as x.
Raises ValueError if either x or y is not finite, and ZeroDivisionError
if x is finite and y is zero.
"""
# Handle error cases.
if not (math.isfinite(x) and math.isfinite(y)):
raise ValueError("x and y must be finite")
if not y:
raise ZeroDivisionError("Division by zero")
# Reduce to the case x and y finite, y positive, x nonnegative.
sign_x, x = sign(x), abs(x)
sign_y, y = sign(y), abs(y)
# If x is *much* larger than y, x / y will overflow. In that case we do a
# form of long division: first divide with a shifted version of x, then
# shift back and continue with the remainder. In the common case where the
# exponents of x and y are not too far apart, the for loop body is executed
# only once. In extreme cases, it'll be executed two or three times.
exp_diff = math.frexp(x)[1] - math.frexp(y)[1]
x_shifts = list(range(exp_diff - 954, 0, -954)) + [0]
acc = 0
for x_shift in x_shifts:
q, r = _fdiv_inner(math.ldexp(x, -x_shift), y)
x = math.ldexp(r, x_shift)
acc += q << x_shift
# Ensure that the returned remainder has the same sign as the original x.
if x < 0.0:
x += y
acc -= 1
x += 0.0
return acc * sign_x * sign_y, x * sign_x
def _fdiv_inner(x: float, y: float) -> tuple[int, float]:
"""
Divide one float by another, giving integer quotient and float remainder.
Assumes y is finite and positive, x is finite, and x / y does not overflow.
Returns a pair (q, r) such that x == q * y + r (exactly) and |r| < |y|.
"""
acc = 0
while y <= abs(x):
q = round(x / y)
# Could use a single fma in place of the two lines below, if supported:
# x := -fma(q, y, -x)
xlo, xhigh = exact_mult(q, y)
x = x - xhigh - xlo
acc += q
return acc, x