Message407098
> All the rounding has already happened at the point where ldexp is called, and the result of the ldexp call is exact.
Sketch of proof:
[Here](https://github.com/python/cpython/blob/4ebde73b8e416eeb1fd5d2ca3283f7ddb534c5b1/Objects/longobject.c#L3929) we have:
shift = Py_MAX(diff, DBL_MIN_EXP) - DBL_MANT_DIG - 2;
from which (assuming IEEE 754 as usual) shift >= -1076. (DBL_MIN_EXP = -1021, DBL_MANT_DIG = 53)
[Here](https://github.com/python/cpython/blob/4ebde73b8e416eeb1fd5d2ca3283f7ddb534c5b1/Objects/longobject.c#L4008) we round away the last two or three bits of x, after which x is guaranteed to be a multiple of 4:
x->ob_digit[0] = low & ~(2U*mask-1U);
Then after converting the PyLong x to a double dx with exactly the same value [here](https://github.com/python/cpython/blob/4ebde73b8e416eeb1fd5d2ca3283f7ddb534c5b1/Objects/longobject.c#L4020) we make the ldexp call:
result = ldexp(dx, (int)shift);
At this point dx is a multiple of 4 and shift >= -1076, so the result of the ldexp scaling is a multiple of 2**-1074, and in the case of a subnormal result, it's already exactly representable.
For the int/int division possibly not being correctly rounded on x87, see [here](https://github.com/python/cpython/blob/4ebde73b8e416eeb1fd5d2ca3283f7ddb534c5b1/Objects/longobject.c#L3889-L3892).
It won't affect _this_ application, but possibly we should fix this anyway. Though the progression of time is already effectively fixing it for us, as x87 becomes less and less relevant. |
|
Date |
User |
Action |
Args |
2021-11-26 21:46:19 | mark.dickinson | set | recipients:
+ mark.dickinson, tim.peters, rhettinger, steven.daprano |
2021-11-26 21:46:19 | mark.dickinson | set | messageid: <1637963179.57.0.0482219202765.issue45876@roundup.psfhosted.org> |
2021-11-26 21:46:19 | mark.dickinson | link | issue45876 messages |
2021-11-26 21:46:19 | mark.dickinson | create | |
|