True division of integers might be invalid due to double rounding

[skirpichev]

Bug report

Bug description:

The long_true_divide() has a fast path for small integers:

cpython/Objects/longobject.c

Lines 4720 to 4740 in 726e8e8

	/* Fast path for a and b small (exactly representable in a double).
	Relies on floating-point division being correctly rounded; results
	may be subject to double rounding on x86 machines that operate with
	the x87 FPU set to 64-bit precision. */
	a_is_small = a_size <= MANT_DIG_DIGITS ||
	(a_size == MANT_DIG_DIGITS+1 &&
	a->long_value.ob_digit[MANT_DIG_DIGITS] >> MANT_DIG_BITS == 0);
	b_is_small = b_size <= MANT_DIG_DIGITS ||
	(b_size == MANT_DIG_DIGITS+1 &&
	b->long_value.ob_digit[MANT_DIG_DIGITS] >> MANT_DIG_BITS == 0);
	if (a_is_small && b_is_small) {
	double da, db;
	da = a->long_value.ob_digit[--a_size];
	while (a_size > 0)
	da = da * PyLong_BASE + a->long_value.ob_digit[--a_size];
	db = b->long_value.ob_digit[--b_size];
	while (b_size > 0)
	db = db * PyLong_BASE + b->long_value.ob_digit[--b_size];
	result = da / db;
	goto success;
	}

As mentioned in the comment, this code relies on correct rounding of floating-point division. Wrong behavior could be e.g. triggered by building the 32-bit CPython on 64-bit Linux host (with CPython floats being IEEE doubles, little-endian):

./configure CFLAGS="-m32" LDFLAGS="-m32"
sed -i -e 's!#define PY_HAVE_PERF_TRAMPOLINE 1!#undef PY_HAVE_PERF_TRAMPOLINE!g' pyconfig.h
make -s

>>> 12143 / 31517  # correct answer is 0.385284132373005
0.3852841323730051

I think it's a bug. Why not enable this part conditionally, if X87_DOUBLE_ROUNDING is not defined? Attached patch fixes the problem for me, it seems that the rest of code doesn't suffer from double rounding and more or less closely follows to the pure-Python version of this routine.

Alternatively, we can even drop special handling for small ints (X87_DOUBLE_ROUNDING may not catch all cases of broken division). I didn't measured yet performance impact of this, but IMO correctness should be more important.

A patch with ifdef X87_DOUBLE_ROUNDING

diff --git a/Objects/longobject.c b/Objects/longobject.c
index 43c0db753a0..dd06e5878dc 100644
--- a/Objects/longobject.c
+++ b/Objects/longobject.c
@@ -4611,7 +4611,7 @@ long_true_divide(PyObject *v, PyObject *w)
     PyLongObject *a, *b, *x;
     Py_ssize_t a_size, b_size, shift, extra_bits, diff, x_size, x_bits;
     digit mask, low;
-    int inexact, negate, a_is_small, b_is_small;
+    int inexact, negate;
     double dx, result;
 
     CHECK_BINOP(v, w);
@@ -4717,14 +4717,15 @@ long_true_divide(PyObject *v, PyObject *w)
     if (a_size == 0)
         goto underflow_or_zero;
 
+#ifndef X87_DOUBLE_ROUNDING
     /* Fast path for a and b small (exactly representable in a double).
        Relies on floating-point division being correctly rounded; results
        may be subject to double rounding on x86 machines that operate with
        the x87 FPU set to 64-bit precision. */
-    a_is_small = a_size <= MANT_DIG_DIGITS ||
+    int a_is_small = a_size <= MANT_DIG_DIGITS ||
         (a_size == MANT_DIG_DIGITS+1 &&
          a->long_value.ob_digit[MANT_DIG_DIGITS] >> MANT_DIG_BITS == 0);
-    b_is_small = b_size <= MANT_DIG_DIGITS ||
+    int b_is_small = b_size <= MANT_DIG_DIGITS ||
         (b_size == MANT_DIG_DIGITS+1 &&
          b->long_value.ob_digit[MANT_DIG_DIGITS] >> MANT_DIG_BITS == 0);
     if (a_is_small && b_is_small) {
@@ -4738,6 +4739,7 @@ long_true_divide(PyObject *v, PyObject *w)
         result = da / db;
         goto success;
     }
+#endif

CPython versions tested on:

CPython main branch

Operating systems tested on:

Linux

[serhiy-storchaka]

Why is it happened in the 32-bit CPython on 64-bit Linux host?

Can the optimization be safely used for smaller integers even in case of double rounding?

[skirpichev]

On Tue, Dec 09, 2025 at 03:16:58AM -0800, Serhiy Storchaka wrote: Why is it happened in the 32-bit CPython on 64-bit Linux host?

Same defaults for FPU precision as for 32-bit box, I think.

Can the optimization be safely used for smaller integers even in case of double rounding?

Well, in case of my system the problem could be also solved by using _Py_SET_53BIT_PRECISION* macros (as in pystrtod.c). But if we don't know how to reset precision on machine, I don't think that this algorithm could be used.

[serhiy-storchaka]

You solution looks correct.

I just wonder if we can enable optimization for really small numbers. Is 31517 the smallest denominator for which you can reproduce the issue? Can we leave the optimization for smaller integers, and can we be sure that it always works?

cc @mdickinson, @tim-one

[mdickinson]

Is 31517 the smallest denominator for which you can reproduce the issue?

There are smaller examples, like 1/2731 or 115/2051. I suspect that it can be proved that a denominator of 2048 or smaller is safe.

It's been a very long time since I've had to work with a machine that does x87-style long double arithmetic, though.

[mdickinson]

Why is it happened in the 32-bit CPython on 64-bit Linux host?

There are two factors at play:

• the 32-bit build likely uses x87 FPU instructions instead of SSE2 instructions
• on Linux, the default x87 precision is typically 64 bits (i.e., long double precision) rather than 53 bits. On Windows I believe it's 53 bits, so even if x87 precision is used you're safe from double rounding, at least for normal numbers. (For subnormals it's another story.)

[mdickinson]

I suspect that it can be proved that a denominator of 2048 or smaller is safe.

Second thoughts: we'd probably definitely also need a bound on the numerator for that to be true. (With unbounded numerators, we can even get double rounding for a denominator of 1.)

[mdickinson]

I think it's a bug.

I'd disagree with this assessment, though. There's currently no promise anywhere that I know of that the result of integer division will be correctly rounded. It's a nice feature to have, for sure.

[tim-one]

It's rare that I disagree with @mdickinson, and this is no exception ;-) It's not "a bug", but would be nice to repair, provided so at reasonable cost. Materially slowing true division of smalls probably wouldn't be "reasonable" to me, although slowing true division of giant ints wouldn't bother me.

I don't have access to the systems at issue, but am confused on one point:

>>> 12143 / 31517  # correct answer is 0.385284132373005
0.3852841323730051

What happens if you do 12143.0 / 31517 0 instead? It's far more important (IMO) that true division of ints return exactly the same result as float division, on a single box, when the ints are exactly representable as floats. Best I understand the issues here, the "all-float" spelling would return the same thing.

[skirpichev]

There are smaller examples, like 1/2731 or 115/2051.

This was triggered by fuzzy testing of float(mpq) (that means mpfr_set_q() with MPFR_RNDN) in the gmpy2 on i686, so example is rather random.

There's currently no promise anywhere that I know of that the result of integer division will be correctly rounded.

Strictly speaking, yes. And I suspect, it's intentional.

On another hand, someone could expect that all methods on integers return exact answers (well defined, independent from hardware).

It's not "a bug", but would be nice to repair, provided so at reasonable cost.

So, in above sense it rather looks as a bug. But maybe you are right and preserve n/d == float(n)/float(d) is more important. Though, there is some lie, as this equivalence will be broken for big enough integers...

About cost. Here is what happens if the special case removed:

Benchmark	ref	patch
12143/31517	180 ns	418 ns: 2.32x slower

Probably it's not an option, as issue is valid for rather exotic systems nowadays.

What happens if you do 12143.0 / 31517.0 instead?

As both operands fit into floating-point type without loss, same answers. (I.e. 0.3852841323730051 on broken systems.)

But we could apply same corrections (_Py_SET_53BIT_PRECISION*) to true division for integers and floats.

diff --git a/Objects/floatobject.c b/Objects/floatobject.c
index 78006783c6e..1b59240654e 100644
--- a/Objects/floatobject.c
+++ b/Objects/floatobject.c
@@ -593,6 +593,7 @@ static PyObject *
 float_div(PyObject *v, PyObject *w)
 {
     double a,b;
+    _Py_SET_53BIT_PRECISION_HEADER;
     CONVERT_TO_DOUBLE(v, a);
     CONVERT_TO_DOUBLE(w, b);
     if (b == 0.0) {
@@ -600,7 +601,9 @@ float_div(PyObject *v, PyObject *w)
                         "division by zero");
         return NULL;
     }
+    _Py_SET_53BIT_PRECISION_START;
     a = a / b;
+    _Py_SET_53BIT_PRECISION_END;
     return PyFloat_FromDouble(a);
 }
 
diff --git a/Objects/longobject.c b/Objects/longobject.c
index 43c0db753a0..981a52d06c2 100644
--- a/Objects/longobject.c
+++ b/Objects/longobject.c
@@ -4613,6 +4613,7 @@ long_true_divide(PyObject *v, PyObject *w)
     digit mask, low;
     int inexact, negate, a_is_small, b_is_small;
     double dx, result;
+    _Py_SET_53BIT_PRECISION_HEADER;
 
     CHECK_BINOP(v, w);
     a = (PyLongObject *)v;
@@ -4735,7 +4736,9 @@ long_true_divide(PyObject *v, PyObject *w)
         db = b->long_value.ob_digit[--b_size];
         while (b_size > 0)
             db = db * PyLong_BASE + b->long_value.ob_digit[--b_size];
+        _Py_SET_53BIT_PRECISION_START;
         result = da / db;
+        _Py_SET_53BIT_PRECISION_END;
         goto success;
     }
 

Edit: Same benchmark for above patch:

Benchmark	ref	patch
12143/31517	180 ns	191 ns: 1.06x slower

[tim-one]

About the docs, they're generally careful not to promise much of anything about how floats behave. When correctly rounded is the promise, they say so (perhaps limited to the decimal module,, and for binary float <-> str conversion, the round() function, some scattered others).

But maybe you are right and preserve n/d == float(n)/float(d) is more important. Though, there is some lie, as this equivalence will be broken for big enough integers...

Except you skipped my crucial qualifier: WHEN n and d are exactly representable as floats. When that's true, there's no coherent mental model the user can construct for "why" n/d and float(n)/float(d) could return different results. But if one or both aren't representable as floats, "oh - just floating-point slop again" will satisfy almost everyone.

I was surprised (and pleased) to see the heroic effort @mdickinson gave to making n/d correctly rounded for huge ints, but as his comments record, it is indeed vulnerable to double rounding when ints are small. In fact, I was surprised at first that giant int / giant int didn't blow up with an overflow error from the obvious float(n)/float(d) implementation.

Anyway, I don't think it's worth much (if any) effort to chase this one. "Extended double" has fallen out of fashion, and HW increasingly doesn't even bother to offer it. In decades, it seems no user has ever even noticed. As an old-school numerical analyst of sorts myself, while I appreciate correct rounding, I really can't get excited about worst-case errors a tiny fraction of a ULP anyway.

[skirpichev]

Except you skipped my crucial qualifier: WHEN n and d are exactly representable as floats.

Not really, this fits to my "big enough" remark ;-)

When that's true, there's no coherent mental model the user can construct for "why" n/d and float(n)/float(d) could return different results.

One such was mentioned: methods on integers are exact, while for floats they might suffer from some hardware quirks. IMO, it's a clear borderline.

But it seems no one willing to push things in this direction. So, I'm going to close this as "won't fix" on the ground:

1. correct rounding is not a documented feature for true division of integers
2. we can't provide a priceless workaround for broken platforms AND keep n/d==float(n)/float(d) identity.