From mboxrd@z Thu Jan  1 00:00:00 1970
Return-Path: <gcc-bugzilla@gcc.gnu.org>
Received: by sourceware.org (Postfix, from userid 48)
	id 993CA3858D37; Fri,  3 Mar 2023 14:28:05 +0000 (GMT)
DKIM-Filter: OpenDKIM Filter v2.11.0 sourceware.org 993CA3858D37
DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gcc.gnu.org;
	s=default; t=1677853685;
	bh=MW24vaGMWpDNZdUBwxEvF4PrdfsRP39tE9MBy2wFf0Q=;
	h=From:To:Subject:Date:In-Reply-To:References:From;
	b=X8l9HK38zBzL0v7vZYil4LmFbRClat+GCLBbwnUtMjDdCoT529DnCY+XBnC38pGFV
	 2DuPHhZCwcvlFQWsx+9TpZl8a8d0EoMuitM7cnMTseu37PTANyUAk1eTIYsbGZMkh4
	 ZDqHVUwSAwQBC9V4bFCe++CXNft3nq4MG4TdTaeA=
From: "jakub at gcc dot gnu.org" <gcc-bugzilla@gcc.gnu.org>
To: gcc-bugs@gcc.gnu.org
Subject: [Bug tree-optimization/109008] [13 Regression] Wrong code in scipy
 package since r13-3926-gd4c2f1d376da6f
Date: Fri, 03 Mar 2023 14:28:05 +0000
X-Bugzilla-Reason: CC
X-Bugzilla-Type: changed
X-Bugzilla-Watch-Reason: None
X-Bugzilla-Product: gcc
X-Bugzilla-Component: tree-optimization
X-Bugzilla-Version: 13.0
X-Bugzilla-Keywords: wrong-code
X-Bugzilla-Severity: normal
X-Bugzilla-Who: jakub at gcc dot gnu.org
X-Bugzilla-Status: NEW
X-Bugzilla-Resolution: 
X-Bugzilla-Priority: P1
X-Bugzilla-Assigned-To: unassigned at gcc dot gnu.org
X-Bugzilla-Target-Milestone: 13.0
X-Bugzilla-Flags: 
X-Bugzilla-Changed-Fields: 
Message-ID: <bug-109008-4-sqlyxUiNtJ@http.gcc.gnu.org/bugzilla/>
In-Reply-To: <bug-109008-4@http.gcc.gnu.org/bugzilla/>
References: <bug-109008-4@http.gcc.gnu.org/bugzilla/>
Content-Type: text/plain; charset="UTF-8"
Content-Transfer-Encoding: quoted-printable
X-Bugzilla-URL: http://gcc.gnu.org/bugzilla/
Auto-Submitted: auto-generated
MIME-Version: 1.0
List-Id: <gcc-bugs.sourceware.org>

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=3D109008

--- 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, doub=
le
*result_max)
{
  *result_min =3D result;
  *result_max =3D result;
  for (double eps =3D __DBL_DENORM_MIN__; __builtin_isfinite (eps); eps =3D=
 eps *
2.0)
    if (!test (result - eps))
      {
        *result_min =3D result - eps;
        break;
      }
  for (double eps =3D __DBL_DENORM_MIN__; __builtin_isfinite (eps); eps =3D=
 eps *
2.0)
    if (!test (result + eps))
      {
        *result_max =3D result + eps;
        break;
      }
}

__attribute__((noipa)) double
fn1 (double eps)
{
  double d =3D 1. + eps;
  if (d =3D=3D 1.)
    return eps;
  return 0.0;
}

int
test1 (double eps)
{
  return fn1 (eps) =3D=3D eps;
}

__attribute__((noipa)) double
fn2 (double eps)
{
  double d =3D 1. + eps;
  if (d =3D=3D 0x1.0000001p0)
    return eps;
  return 0.0;
}

int
test2 (double eps)
{
  return fn2 (eps) =3D=3D eps;
}

__attribute__((noipa)) double
fn3 (double eps)
{
  double d =3D 1. + eps;
  if (d =3D=3D 2.)
    return eps;
  return 0.0;
}

int
test3 (double eps)
{
  return fn3 (eps) =3D=3D 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. =3D=3D 2. reverse is affected, but while ma=
ximum of
other_op and lhs range's ulp might be a good starting point, it would be go=
od
to iteratively adjust it.  The above just searches for a power of two epsil=
on
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 bi=
nary
search between the largest and smallest exponents,
and/or whether after finding a pair of power of two epsilons where one exte=
nds
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 aft=
er
some constant number of them and consider that good enough (similarly, if t=
he
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 a=
nd
do the math right.=