[Bug tree-optimization/126562] New: Invalid pow (x, 0.5) -> sqrt transformation when signed zeros or infinities are required

[ktkachov at gcc dot gnu.org via Gcc-bugs]

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

            Bug ID: 126562
           Summary: Invalid pow (x, 0.5) -> sqrt transformation when
                    signed zeros or infinities are required
           Product: gcc
           Version: 17.0
            Status: UNCONFIRMED
          Keywords: wrong-code
          Severity: normal
          Priority: P3
         Component: tree-optimization
          Assignee: unassigned at gcc dot gnu.org
          Reporter: ktkachov at gcc dot gnu.org
  Target Milestone: ---

I guess these are two separate bugs but they are related:

/* pow (x, 0.5) must not become sqrt (x) while signed zeros are honoured:
   C99 Annex F.10.4.4 gives pow (-0.0, 0.5) == +0.0, while F.10.4.5 gives
   sqrt (-0.0) == -0.0.

   Two independent sites did this rewrite with an insufficient guard:
     tree-vect-patterns.cc vect_recog_pow_pattern, no guard at all
     tree-ssa-math-opts.cc gimple_expand_builtin_pow, no !HONOR_INFINITIES */

#include <math.h>

void __attribute__((noinline, noipa))
f (double *restrict d, const double *restrict s, int n)
{
  for (int i = 0; i < n; i++)
    d[i] = pow (s[i], 0.5);
}

int
main (void)
{
  double s[8], d[8];

  for (int i = 0; i < 8; i++)
    s[i] = -0.0;

  f (d, s, 8);

  for (int i = 0; i < 8; i++)
    if (__builtin_signbit (d[i]) != 0)
      __builtin_abort ();

  return 0;
}

/* pow (x, 0.5) must not become sqrt (x) while infinities are honoured:
   C99 Annex F.10.4.4 gives pow (-Inf, 0.5) == +Inf, while F.10.4.5 gives
   sqrt (-Inf) == NaN, with the invalid exception raised.

   The comment on the scalar fold in tree-ssa-math-opts.cc claims the
   rewrite is "always safe unless signed zeros must be maintained", which
   misses this case entirely. */

#include <math.h>

void __attribute__((noinline, noipa))
f (double *restrict d, const double *restrict s, int n)
{
  for (int i = 0; i < n; i++)
    d[i] = pow (s[i], 0.5);
}

int
main (void)
{
  double s[8], d[8];

  for (int i = 0; i < 8; i++)
    s[i] = 4.0;
  s[3] = -__builtin_inf ();

  f (d, s, 8);

  if (d[3] != __builtin_inf ())
    __builtin_abort ();

  return 0;
}

Both abort at -O3 -fno-math-errno on aarch64