[Bug fortran/96811] New: Power: x**(negative integer) – use libgcc variant for power raised to negative value?

[burnus at gcc dot gnu.org]

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

            Bug ID: 96811
           Summary: Power: x**(negative integer) – use libgcc variant for
                    power raised to negative value?
           Product: gcc
           Version: 11.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: fortran
          Assignee: unassigned at gcc dot gnu.org
          Reporter: burnus at gcc dot gnu.org
  Target Milestone: ---

Fabio Azevedo has some observations at
https://gcc.gnu.org/pipermail/fortran/2020-August/054931.html

In particular:

  The only difference in behaviour I see is that when the exponent is
  negative, libgcc2 calculates 1/(x**(-n)) while libgfortran calculates
  (1/x)**(-n). The first version being more accurate.
[...]
  On the top of that x**n with x real for negative exponents appears to yield
  results more compatible with 1/(x**(-n)) than (1/x)**(-n).

libgcc/libgcc2.c has:

NAME (TYPE x, int m)
{
  unsigned int n = m < 0 ? -m : m;
  TYPE y = n % 2 ? x : 1;
  while (n >>= 1)
    {
      x = x * x;
      if (n % 2)
        y = y * x;
    }
  return m < 0 ? 1/y : y;  // <<< division done at the end.
}


While libgfortran/m4/pow.m4 has (for noninteger x):

      if (n < 0)
        {
          u = -n;
          x = pow / x;  // <<< division done already here; pow == 1
`       }
      else
        {
           u = n;
        }
      for (;;)
        {
          if (u & 1)
            pow *= x;
          u >>= 1;
          if (u)
            x *= x;
          else
            break;
        }
    }
  return pow;