Numerical methods: why doesn't this python code return 1.0?

Question

I typed the following into the python console:

>>>from numpy import float64
>>>x=float64(1.98682855148322934369)
>>>x
np.float64(1.9868285514832293)
>>>y=float64(1)/x
>>>x*y
np.float64(0.9999999999999999)

My argument for why the above code should return 1.0 is as follows:

1. The numbers follow IEEE double precision, meaning that its $\epsilon_{mach}$ (for round to nearest) is $2^{-53}$ (i.e. the minimum $\epsilon$ so that $fl(1+\epsilon)\ne1$)
2. We know 1/x is in the range $(2^{-1},1)$, so that the maximum absolute error for fl(1/x) is $2^{-1}*2^{-53}=2^{-54}$
3. Let $fl(1/x)=1/x+\epsilon$ where $\epsilon$ is the error with absolute value $<=2^{-54}$
4. Then $x*fl(1/x)=1+x\epsilon$. Note that $x<2$, so $|x\epsilon|<2*2^{-54}=2^{-53}=\epsilon_{mach}$
5. Therefore, $fl(x*fl(1/x))=1$

I checked my reasoning with deepseek r1 and it agreed with me. I also performed all sorts of hardware checks to make sure all intermediate steps were 64 bit precision and everything seemed to be fine. So this is perplexing. Why didn't python return 1 as the answer?

Answer 1

First, the float64 stuff isn't necessary - Python floats are already float64s.

You've got correctly that the difference $x\epsilon$ between $xy$ and 1 is at most $2^{-53}$, but that's not enough to say that the nearest float to $xy$ is 1. Here $\epsilon$ will be negative (it's -0x1.1956...p-55) and the float immediately before 1 is $1-2^{-53}$, so the halfway point between 1 and the preceding float is $1-2^{-54}$, and just knowing $|x\epsilon | \leqslant 2^{-53}$ isn't enough to say which side of this halfway point $1+x\epsilon$ is. In fact $x\epsilon$ is -0x1.177...p-54, larger in magnitude than $2^{-54}$, so the nearest float to $xy$ will be $1-2^{-53}$.

The mpmath arbitrary precision library is useful for experimenting with this sort of thing, as is numpy.nextafter(a, b) which returns the next float after a in the direction of b.

Answer 2

5 does not follow from 4. The error must be less than half of the resolution for 1 to be the closest representable number.