BUG: Fix <complex 0>^{non-zero} by prithvitewatia · Pull Request #18535 · numpy/numpy

prithvitewatia · 2021-03-03T13:57:36Z

The IEEE standard 754 2008 is used in the development of floating-point libraries.
The 0^0 case is discussed in Wikipedia.
For the case of 0^non-zero, the result is the same as that of C99.
Side note: C99 for negative real in complex powers generates inf + nan but we are generating nan + nan as
it is consistent with our previous warning.
Resolves #18378.

numpy/core/src/npymath/npy_math_complex.c.src

seberg · 2021-03-03T17:34:48Z

numpy/core/tests/test_umath.py

+            #real part less than zero imag part greater than zero
+            assert_complex_equal(np.power(zero,-1+1j),cnan)
+            #real part less than zero imag part less than zero
+            assert_complex_equal(np.power(zero,-2-3j),cnan)


Can you split these out and use pytest.warns()? That way it is explicitly about which case gives warnings and which does not. (and the test will notify it if it does not work properly due to a change or weird platform)

prithvitewatia · 2021-03-06T02:34:47Z

@eric-wieser @seberg I have made the changes. Please review my pull request.

prithvitewatia · 2021-03-10T16:29:38Z

@eric-wieser @seberg friendly reminder. I don't mean to bother you but you haven't reviewed my pull request since the last changes. Please take a look.

seberg · 2021-03-10T16:49:36Z

A couple of notes. Sorry, I have to dig into some things or I won't make progress there. I thin kthis goes in the right direction, but I wonder if we can't simplify the code quite a bit.

Could you try to make sure to follow the C and python style guide. I.e. no missing spaces, multi-line comments styled as:
```
 /* 
  * Thi sis a multiline
  * comment
  */
```
Can you check if assuming #ifdef HAVE_CPOW@C@ is always true is OK? It should be, since cpow is part of C99.
I am surprised that the integer special case is worthwhile for up to -100 to +100! (side-note)
It would be cool if we could skip these "corner cases" in the fast-path if block and just fall through to the system cpow (see point 2.). Then we don't have to worry about it (although chances are that clang will come around and be buggy...)
Don't do warning captures of tests that don't produce warnings, it is unnecessary. Use match= with pytest.warns if you want to check the warning message context. It is typically unnecessary to use with pytest.warns(...) as rec

prithvitewatia · 2021-03-11T09:12:37Z

@seberg I have checked cpow in C99 and C11 and both are present and working fine. Also could you please explain your third point as I am not able to understand it?

seberg · 2021-03-11T10:21:44Z

Cool, thanks, sorry very briefly. This code block optimizes integer powers:

numpy/numpy/core/src/npymath/npy_math_complex.c.src

Lines 470 to 506 in a157224

    
           if (bi == 0 && (n=(npy_intp)br) == br) { 
        
               if (n == 1) { 
        
                   /* unroll: handle inf better */ 
        
                   return npy_cpack@c@(ar, ai); 
        
               } 
        
               else if (n == 2) { 
        
                   /* unroll: handle inf better */ 
        
                   return cmul@c@(a, a); 
        
               } 
        
               else if (n == 3) { 
        
                   /* unroll: handle inf better */ 
        
                   return cmul@c@(a, cmul@c@(a, a)); 
        
               } 
        
               else if (n > -100 && n < 100) { 
        
                   @ctype@ p, aa; 
        
                   npy_intp mask = 1; 
        
                   if (n < 0) { 
        
                       n = -n; 
        
                   } 
        
                   aa = c_1@c@; 
        
                   p = npy_cpack@c@(ar, ai); 
        
                   while (1) { 
        
                       if (n & mask) { 
        
                           aa = cmul@c@(aa,p); 
        
                       } 
        
                       mask <<= 1; 
        
                       if (n < mask || mask <= 0) { 
        
                           break; 
        
                       } 
        
                       p = cmul@c@(p,p); 
        
                   } 
        
                   r = npy_cpack@c@(npy_creal@c@(aa), npy_cimag@c@(aa)); 
        
                   if (br < 0) { 
        
                       r = cdiv@c@(c_1@c@, r); 
        
                   } 
        
                   return r; 
        
               }

And the last part optimizes all between -100 and 100, although I can't say I understand the algorithm, so it probably makes sense :).

prithvitewatia · 2021-03-15T18:16:48Z

@seberg friendly reminder. I have made the changes please take a look.

charris · 2021-03-26T18:26:54Z

close/reopen

eric-wieser · 2021-03-26T18:34:34Z

Can you update the PR description to reference the bit of the IEEE standard that says this is correct behavior?

prithvitewatia · 2021-03-27T12:45:06Z

@eric-wieser I have made the changes. Please take a look.

prithvitewatia · 2021-03-27T14:41:16Z

@seberg friendly reminder. Please review my pull request.

prithvitewatia · 2021-04-04T17:22:09Z

@seberg I apologize for tagging you again and again but you haven't reviewed my pull request. I just wanted to remind you in case it slipped.

prithvitewatia · 2021-04-17T18:10:47Z

@seberg friendly reminder. Please review my Pull request.

seberg

Sorry, I am still confused about the warning flags a bit. And now also about the +0 and -0 dependency, which the (deleted) NB: comment allures too.

numpy/core/src/npymath/npy_math_complex.c.src

seberg · 2021-04-28T17:37:33Z

numpy/core/src/npymath/npy_math_complex.c.src

+        return npy_cpack@c@(1.,0.);
+    }
+    /*Here it is already tested that for a^b, b is not zero.
+    So we are checking behaviour of a.*/


I don't think thois comment is necessary? It might make sense in the "generating an invalid" branch below.

@seberg Could you please explain what do you mean by generating an invalid as I am not able to understand it

I mean moving it down to the branch where you return NaN manually. Maybe its a bit of a code layout thing in the diff, but I think that branch is the place where you may want to point out that 0^0 is already handled.

All right got it. Thanks

numpy/core/src/npymath/npy_math_complex.c.src

seberg · 2021-04-28T17:43:50Z

numpy/core/src/npymath/npy_math_complex.c.src

+            //Generating an invalid
+            volatile @type@ tmp=NPY_INFINITY@C@;
+            tmp-=NPY_INFINITY@C@;
+            ar=tmp;


Please fix the spaces here again, no operators without spaces around them.

@seberg Here is what I found.
The first section discusses how the complex exponentiation is done.
Second section deals with the sign.

Complex Exponentiation
If z and w are two complex numbers then
z^w = e^{w*(ln|z|+i*arg(z} *e^{i*w*(2*n*pi)}
which means that z^w is multivalued. We call the value at n=0 the principle value.
Now, here I will discuss the case of 0^{non-zero complex number} which is related to our issue.
0^w = e^{w * (ln(0) + i*arg(0,0))}*e^{ i*w*(2*n *pi)}
Lets call this equation 1.
Now arg is calculated with the help of atan2. Mathematically undefined but is realized in computer languages.
As we are considering C language, for unsigned 0, we get the output as 0. For signed numbers, we get the
output 0, pi, -0, -pi. This is demonstrated with the following C program.

#include<stdio.h> #include<math.h> int main(){ double a=atan2(0,0); double b=atan2(+0.0,+0.0); double c=atan2(+0.0,-0.0); double d=atan2(-0.0,+0.0); double e=atan2(-0.0,-0.0); printf("%lf\n%lf\n%lf\n%lf\n%lf\n",a,b,c,d,e); return 0; }

Output

0.000000 0.000000 3.141593 -0.000000 -3.141593

Therefore e^{i*arg(0,0)} will be +1 or -1 depending upon values of ar and ai (0 but still signed).
Now, let's spilt equation 1 into three parts.
Part 1: e^{ w*ln|0| } (ln|0| = -inf)
Part 2: e^{ i*arg(0,0) }= ±1
Part 3: e^{ i*w*(2n*pi)}
For complex number w , we get 4 cases
br = 0, bi != 0
br != 0, bi = 0
br = 0, bi != 0
br != 0, bi != 0
The case for br=0 and bi=0 leads to 0^0 which is defined to be 1 by IEEE.
For case br != 0, bi = 0, we have

if br > 0 we have Part 1= e^{ br* -inf }= 0 Part2 = ±1 (actual sign depending on signs of ar and ai) Part3 = dont care Result = ±0 if br < 0 we have Part 1= e^{ br*-inf }= inf Part 2= ±1 Part 3= oscillates between positive and negative value Therefore this leads to indeterminate form of inf*0 and the result is undefined . (This might also expalin why C generates (inf,nan) as its output when encountered with such values. In our case we generate (nan,nan).

For br = 0, bi != 0

if bi > 0 Part 1: e^{ i*bi*-inf} which can be written as cos(bi*-inf)+isin(bi*-inf) which is undefined. Part 2: don't care Part 3: don't care Result: undefined if bi < 0 Same as above because this time we get +inf instead of minus inf.

For br !=0 and bi != 0

Part 1: e^{ w * -inf} =e^{br*-inf} * e^{i*bi*-inf} (w=br + i*bi) Therefore if br > 0 this goes down to 0 if not, the result oscillates and is undifined. Part 2 and 3 are dont care.

All the cases are summarized in the table below

bi < 0 bi = 0 bi > 0

br < 0 undef undef undef

br = 0 undef 1 undef

br > 0 0 0 0

Value of 0 ^ {complex}

Note: 0s are signed.

Combining with the complex exponentiation table

exp = 0 exp = non-zero

mantissa = 0 1 Table above

mantissa = non-zero 1 Already handled

As we learned that complex exponentiation is multi-valued and the actual
signs depend on many factors such as signs of ar and ai and whether the API
is selecting principle value or not. Also 0+0i, 0-0i,-0+0i and -0-0i correspond
to the same point on the complex plane. So as long the result is 0 it mathematically
does not affect the result

I am ready to go with whatever you suggest. If you suggest going with signs (the resulting 0 contains sign opposite of bi as the C program below shows.)
I will make changes in the code accordingly. If you say its fine just make the code according to numpy C style I will do that.

#include<stdio.h> #include<math.h> #include<complex.h> int main(){ complex mantissa[]={ 0.0+0.0*I, 0.0-0.0*I, -0.0+0.0*I, -0.0-0.0*I}; 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++){ for(int j=0;j<4;j++){ complex c=cpow(mantissa[i],exponents[j]); printf("%lf + %lf *I\n",creal(c),cimag(c)); } } return 0; }

Output:

0.000000 + -0.000000 *I 0.000000 + 0.000000 *I inf + -nan *I inf + -nan *I 0.000000 + -0.000000 *I 0.000000 + 0.000000 *I inf + -nan *I inf + -nan *I 0.000000 + -0.000000 *I 0.000000 + 0.000000 *I inf + -nan *I inf + -nan *I 0.000000 + -0.000000 *I 0.000000 + 0.000000 *I inf + -nan *I inf + -nan *I

References:
Complex Eponentiation
atan2

Sorry for taking a while to respond. That anlysis sounds great, and atan2 is a seems like a great "basis" ~~to base~~ for a decision, I think, thanks! (I have not followed every step, I admit)

My feeling is that we should try to preserve the signs for the 0 return. And the fact that GCC gives a floating point warnings seems like a small issue in glibc? After all, atan2 makes a point of getting this right as well. Also the old code-comment to me reads like: Let's return NaN, because the sign means that returning 0 would not be quite right.

So my personal judgment would be: return the correct sign here, in the special case! And IIRC the rest already looked good. I double checked the C99 standard, but could not find anything relevant (the cpow documentation is pretty much a stub).

I think we should open an issue on gcc/glibc about the warning flags being set for those signed 0 returns. In the best case scenario we learn something new. In the other best case scenario glibc stops giving spurious warnings ;).

Ah GLIBC issue of course, I have opened a bug report, hopefully in the right place here: https://sourceware.org/bugzilla/show_bug.cgi?id=27875

seberg · 2021-04-28T17:45:10Z

Otherwise, the main thing is that it would be nice if you can clean up the code to adher to NumPy C-style.

prithvitewatia · 2021-05-24T07:24:58Z

@seberg Friendly reminder. I have made the changes. Please review my Pull request.

seberg

Just some nits and the tightening the tests. Otherwise, I think this is good to go in. Thanks!

numpy/core/src/npymath/npy_math_complex.c.src

numpy/core/tests/test_umath.py

numpy/core/src/npymath/npy_math_complex.c.src

seberg

Thanks, looks good to me now. I am considering if we should add a release note, but in general, I think it can go in as-is.

prithvitewatia · 2021-06-09T07:03:52Z

@seberg I think we should add a release note (in the changes section or improvements section whichever you suggest). Please guide me on how to add a release note and I will add it. Many thanks. 😃

seberg · 2021-06-09T13:48:31Z

You have to add a brief rst file to doc/release/upcoming_changes. The README.txt lists the categories. Probably be an improvement or compatibility if you think users might actually notice. (There is currently only one example because 1.21 was just branched.)

seberg · 2021-06-16T20:03:19Z

I had meant to just fix up the release notes myself a bit (pushed a change, feel free to overwrite).

But now I noticed that we are not handling -0. on the inputs :(. And I am not quite sure what is right. Also, the quality check in the tests do not actually check for -0., we would have to use np.signbit there to check for the -0. output.

prithvitewatia · 2021-06-18T02:20:16Z

@seberg Thanks for enhancing the release note 😄 . I also thought of adding this note but then thought this is documentation related so I didn't add added it in the release note.
I am sorry I also missed that we are not checking at negative zero in the inputs. I think that all the tests should be done for four cases of zero input in the mantissa { (0., 0.), (0.,-0.),(-0.,0.),(-0.,-0.) }.
I will fix these issues in the next commit.
Also, I will check for -0. using np.signbit.
Many thanks for reminding me 🤝 .

Co-authored-by: Ivan Gonzalez <scratchmex@gmail.com>

Co-authored-by: Matti Picus <matti.picus@gmail.com> Co-authored-by: Ivan Gonzalez <scratchmex@gmail.com>

prithvitewatia · 2022-06-17T15:47:49Z

The circleci problem should go away when merged or rebased. @prithvitewatia gentle ping.

Hi @charris, sorry for the delay. I've committed the suggestions and rebased with the latest main now.

charris · 2022-11-19T18:23:19Z

@seberg Are you good with this?

seberg · 2022-11-21T11:23:22Z

Sorry, I am burned on this. I would like someone else to make up their mind carefully and then, yes I am happy. But I am not Kahan and all standards I am aware of bail on this beyond saying cexp(c*clog(x)) is an acceptable implementation (meaning that e.g. warnings are allowed to be off in these corner cases).
Right now, the status is that we do not return 0-0j (signed zeros) unlike a previous version. The release note seems to still refer to the previous version.

charris · 2022-11-21T13:27:14Z

I am burned on this.

Sorry :) I just wanted to get it off the 1.24 checklist, and away it goes.

doc/release/upcoming_changes/18535.improvement.rst

numpy/core/src/npymath/npy_math_complex.c.src

mattip · 2022-11-30T16:55:04Z

I punted on all the places this says the value is defined by IEEE 754 since we have not seen such a definition, and think we should put this in since it replaces an unjustified NAN result.

Co-authored-by: Matti Picus <matti.picus@gmail.com>

prithvitewatia · 2022-12-01T17:12:14Z

@mattip gentle ping. I have updated the PR according to changes suggested. Please take a look.

mattip · 2022-12-02T11:41:23Z

Thanks @prithvitewatia

github-actions bot added the 00 - Bug label Mar 3, 2021

eric-wieser reviewed Mar 3, 2021

View reviewed changes

numpy/core/src/npymath/npy_math_complex.c.src Show resolved Hide resolved

seberg reviewed Mar 3, 2021

View reviewed changes

Base automatically changed from master to main March 4, 2021 02:05

seberg self-requested a review March 15, 2021 18:47

charris closed this Mar 26, 2021

charris reopened this Mar 26, 2021

seberg reviewed Apr 28, 2021

View reviewed changes

prithvitewatia requested a review from seberg May 15, 2021 03:39

seberg reviewed May 24, 2021

View reviewed changes

seberg self-requested a review June 2, 2021 21:17

seberg approved these changes Jun 7, 2021

View reviewed changes

prithvitewatia and others added 10 commits June 17, 2022 21:13

Fixed Tests with pytest.warns

1587097

Fixed style and warning capture

7c80b55

Fixed whitespace after comma in test_umath.py

9c91f21

Fixed 0 sign and added new tests

890ad66

Fixed spacing and check num of warnings in tests

bf432e2

Added improvement release note

f75ad5e

Fixup release notes a bit

50b3e44

add 754 to IEEE 754

54cbbfa

Co-authored-by: Ivan Gonzalez <scratchmex@gmail.com>

rephrase 0^b comment

4f0d592

Co-authored-by: Ivan Gonzalez <scratchmex@gmail.com>

Apply suggestions from code review

2805164

Co-authored-by: Matti Picus <matti.picus@gmail.com> Co-authored-by: Ivan Gonzalez <scratchmex@gmail.com>

prithvitewatia force-pushed the Issue18378 branch from f067731 to 2805164 Compare June 17, 2022 15:44

InessaPawson changed the title ~~BUG: Fixed <complex 0>^{non-zero} as discussed in issue #18378~~ BUG: Fix <complex 0>^{non-zero} Nov 2, 2022

charris modified the milestones: 1.24.0 release, 1.25.0 release Nov 21, 2022

mattip reviewed Nov 30, 2022

View reviewed changes

doc/release/upcoming_changes/18535.improvement.rst Outdated Show resolved Hide resolved

mattip reviewed Nov 30, 2022

View reviewed changes

doc/release/upcoming_changes/18535.improvement.rst Outdated Show resolved Hide resolved

mattip reviewed Nov 30, 2022

View reviewed changes

numpy/core/src/npymath/npy_math_complex.c.src Outdated Show resolved Hide resolved

prithvitewatia and others added 3 commits December 1, 2022 22:28

Update doc/release/upcoming_changes/18535.improvement.rst

cc919b5

Co-authored-by: Matti Picus <matti.picus@gmail.com>

Update numpy/core/src/npymath/npy_math_complex.c.src

4c55e1b

Co-authored-by: Matti Picus <matti.picus@gmail.com>

Update doc/release/upcoming_changes/18535.improvement.rst

1b57a63

Co-authored-by: Matti Picus <matti.picus@gmail.com>

mattip merged commit 9989e1e into numpy:main Dec 2, 2022

mattip mentioned this pull request Dec 2, 2022

BLD: enscriptem CI run failed #22712

Closed

seiko2plus mentioned this pull request Dec 19, 2022

BUG, SIMD: Fix invalid value encountered in several ufuncs #22834

Merged

Uh oh!

Conversation

prithvitewatia commented Mar 3, 2021 • edited by InessaPawson Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

prithvitewatia commented Mar 6, 2021

Uh oh!

prithvitewatia commented Mar 10, 2021

Uh oh!

seberg commented Mar 10, 2021

Uh oh!

prithvitewatia commented Mar 11, 2021

Uh oh!

seberg commented Mar 11, 2021

Uh oh!

prithvitewatia commented Mar 15, 2021

Uh oh!

charris commented Mar 26, 2021

Uh oh!

eric-wieser commented Mar 26, 2021

Uh oh!

prithvitewatia commented Mar 27, 2021

Uh oh!

prithvitewatia commented Mar 27, 2021

Uh oh!

prithvitewatia commented Apr 4, 2021

Uh oh!

prithvitewatia commented Apr 17, 2021

Uh oh!

seberg left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

seberg May 17, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

seberg commented Apr 28, 2021

Uh oh!

prithvitewatia commented May 24, 2021

Uh oh!

seberg left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

seberg left a comment

Choose a reason for hiding this comment

Uh oh!

prithvitewatia commented Jun 9, 2021

Uh oh!

seberg commented Jun 9, 2021

Uh oh!

seberg commented Jun 16, 2021

Uh oh!

prithvitewatia commented Jun 18, 2021

Uh oh!

prithvitewatia commented Jun 17, 2022

Uh oh!

prithvitewatia commented Mar 3, 2021 •

edited by InessaPawson

Loading

seberg May 17, 2021 •

edited

Loading