https://gcc.gnu.org/bugzilla/show_bug.cgi?id=109008
--- Comment #12 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
(In reply to Richard Biener from comment #11)
I think before we code something on the compiler side, it might be better
to play just with C testcases.
Quite naive
#include <math.h>
__attribute__((noipa)) void
extend_range (double result, int (*test) (double), double *result_min, double
*result_max)
{
*result_min = result;
*result_max = result;
for (double eps = __DBL_DENORM_MIN__; __builtin_isfinite (eps); eps = eps *
2.0)
if (!test (result - eps))
{
*result_min = result - eps;
break;
}
for (double eps = __DBL_DENORM_MIN__; __builtin_isfinite (eps); eps = eps *
2.0)
if (!test (result + eps))
{
*result_max = result + eps;
break;
}
}
__attribute__((noipa)) double
fn1 (double eps)
{
double d = 1. + eps;
if (d == 1.)
return eps;
return 0.0;
}
int
test1 (double eps)
{
return fn1 (eps) == eps;
}
__attribute__((noipa)) double
fn2 (double eps)
{
double d = 1. + eps;
if (d == 0x1.0000001p0)
return eps;
return 0.0;
}
int
test2 (double eps)
{
return fn2 (eps) == eps;
}
__attribute__((noipa)) double
fn3 (double eps)
{
double d = 1. + eps;
if (d == 2.)
return eps;
return 0.0;
}
int
test3 (double eps)
{
return fn3 (eps) == eps;
}
int
main ()
{
double min, max;
extend_range (1. - 1., test1, &min, &max);
__builtin_printf ("%.32a [%.32a, %.32a]\n", 1.0 - 1., min, max);
extend_range (0x1.0000001p0 - 1., test2, &min, &max);
__builtin_printf ("%.32a [%.32a, %.32a]\n", 0x1.0000001p0 - 1., min, max);
extend_range (2. - 1., test3, &min, &max);
__builtin_printf ("%.32a [%.32a, %.32a]\n", 2. - 1., min, max);
return 0;
}
prints
0x0.00000000000000000000000000000000p+0
[-0x1.00000000000000000000000000000000p-53,
0x1.00000000000000000000000000000000p-52]
0x1.00000000000000000000000000000000p-28
[0x1.fffffe00000000000000000000000000p-29,
0x1.00000100000000000000000000000000p-28]
0x1.00000000000000000000000000000000p+0
[0x1.ffffffffffffe0000000000000000000p-1,
0x1.00000000000020000000000000000000p+0]
so yes, clearly even the x + 1. == 2. reverse is affected, but while maximum of
other_op and lhs range's ulp might be a good starting point, it would be good
to iteratively adjust it. The above just searches for a power of two epsilon
in both directions (first one that changes the range, so conservatively
larger), wonder if for the search instead of iterations we couldn't do a binary
search between the largest and smallest exponents,
and/or whether after finding a pair of power of two epsilons where one extends
the lhs range and the other one still doesn't do up to mantissa precision
further iterations to find the exact epsilon value; though perhaps stop after
some constant number of them and consider that good enough (similarly, if the
min - eps (or max + eps) for two different values yield the same value).
All this still sounds quite brute force as opposed to trying to be clever and
do the math right.
