[Bug math/27875] New: Complex power `0^1` gives good result but sets division error flag

[Bug math/27875] New: Complex power `0^1` gives good result but sets division error flag

[sebastian-glibc at sipsolutions dot net]

https://sourceware.org/bugzilla/show_bug.cgi?id=27875

            Bug ID: 27875
           Summary: Complex power `0^1` gives good result but sets
                    division error flag
           Product: glibc
           Version: 2.31
            Status: UNCONFIRMED
          Severity: normal
          Priority: P2
         Component: math
          Assignee: unassigned at sourceware dot org
          Reporter: sebastian-glibc at sipsolutions dot net
  Target Milestone: ---

We had a request in NumPy to return a non NaN value for `0j^(1 + Xj)`. 
Exploring that, I think that glibc gives perfectly good values, including the
correct sign of the 0.

However, the divide by zero floating error flag is set.  Since the return value
is finite, I would not expect it to get set?

This is the result of a test code run (note the first two examples):

```
Calculating for (0 + 0j)^(1 + 1j):
    Result is: 0+-0j
    Divide by zero flag set.

Calculating for (0 + 0j)^(1 + -1j):
    Result is: 0+0j
    Divide by zero flag set.

Calculating for (0 + 0j)^(-1 + 1j):
    Result is: inf+-nanj
    Divide by zero flag set.
    Invalid flag set.

Calculating for (0 + 0j)^(-1 + -1j):
    Result is: inf+-nanj
    Divide by zero flag set.
    Invalid flag set.
```

Generated with the following program:

```
#include <stdio.h>
#include <fenv.h>
#include <complex.h>


static void
print_cpow(complex mantissa, complex exponent)
{
   printf("\nCalculating for (%g + %gj)^(%g + %gj):\n",
          creal(mantissa), cimag(mantissa), creal(exponent), cimag(exponent));
   fetestexcept(FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW | FE_INVALID);

   complex res = cpow(mantissa, exponent);

   int except = fetestexcept(FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW |
FE_INVALID);

   printf("    Result is: %g+%gj\n", creal(res), cimag(res));
   if (except & FE_DIVBYZERO) {
        printf("    Divide by zero flag set.\n");
   }
   if (except & FE_OVERFLOW) {
        printf("    Overflow flag set.\n");
   }
   if (except & FE_UNDERFLOW) {
       printf("    Underflow flag set.\n");
   }
   if (except & FE_INVALID) {
       printf("    Invalid flag set.\n");
   }
}


int
main(void)
{
    complex exponents[] = {1.0+1.0*I, 1.0-1.0*I, -1.0+1.0*I, -1.0-1.0*I};

    for(int i = 0; i < 4; i++){
        print_cpow(-0, exponents[i]);
    }
    return 0;
}
```

And for completeness:

* gcc (Debian 10.2.1-6) 10.2.1 20210110
* ldd (Debian GLIBC 2.31-11) 2.31
* Compiled without any optimization flags: gcc cexp.c -lm

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/27875] Complex power `0^1` gives good result but sets division error flag

[schwab@linux-m68k.org]

https://sourceware.org/bugzilla/show_bug.cgi?id=27875

--- Comment #1 from Andreas Schwab <schwab@linux-m68k.org> ---
cpow is allowed to raise spurious exceptions.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/27875] Complex power `0^1` gives good result but sets division error flag

[schwab@linux-m68k.org]

https://sourceware.org/bugzilla/show_bug.cgi?id=27875

--- Comment #2 from Andreas Schwab <schwab@linux-m68k.org> ---
Specifically, cpow(x,c) raises the same exceptions as cexp(c*clog(x)).  In that
context, a divide-by-zero exception is not spurious.

-- 
You are receiving this mail because:
You are on the CC list for the bug.