On 12/03/19 22:49 +0000, Joseph Myers wrote:
On Tue, 5 Mar 2019, Jonathan Wakely wrote:
The midpoint and lerp functions for floating point types come straight
from the P0811R3 proposal, with no attempt at optimization.
I don't know whether P0811R3 states different requirements from the public
P0811R2, but the implementation of midpoint using isnormal does *not*
satisfy "at most one inexact operation occurs" and is *not* correctly
rounded, contrary to the claims made in P0811R2.
I did wonder how the implementation in the paper was meant to meet the
stated requirements, but I didn't wonder too hard.
Consider e.g. midpoint(DBL_MIN + DBL_TRUE_MIN, DBL_MIN + DBL_TRUE_MIN).
The value DBL_MIN + DBL_TRUE_MIN is normal, but dividing it by 2 is
inexact (and so that midpoint implementation would produce DBL_MIN as
result, so failing to satisfy midpoint(x, x) == x).
Replacing isnormal(x) by something like isgreaterequal(fabs(x), MIN*2)
would avoid those inexact divisions, but there would still be spurious
overflows in non-default rounding modes for e.g. midpoint(DBL_MAX,
DBL_TRUE_MIN) in FE_UPWARD mode, so failing "No overflow occurs" if that's
meant to apply in all rounding modes.
Thanks for this review, and the useful cases to test. Ed is working on
adding some more tests, so maybe he can also look at improving the
code :-)
