bugGNU Octave - Bugs: bug #54302, power operator: 0.^0 results in...

 
 

bug #54302: power operator: 0.^0 results in NaN in some case

Submitted by:  None
Submitted on:  Fri 13 Jul 2018 07:22:49 AM UTC  
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Confirmed Assigned to:  None
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  Open Release:  4.4.0
Operating System:  GNU/Linux

Add a New Comment(Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

( Jump to the original submission)

Mon 16 Jul 2018 03:44:34 AM UTC, comment #6:

I was just thinking back, regarding the

that I removed in the patch, is it possibly the case that at one time these routines used

and that was supposed to find the version of pow() I listed in the previous post, originating from oct-inttypes.cc? That integer version might be slightly more accurate than std::pow(), but possibly a little slower...not drastically slow because it does use a fairly fast reduction by squaring the running product and bit-shifting b. I would imagine that std::pow() is using logarithms or something similar to turn power into multiplication.

Dan Sebald <sebald>
Mon 16 Jul 2018 12:09:03 AM UTC, comment #5:

I have made a modification that seems solve the problem. I basically re-implemented all seven cases of std::pow and pow from the C++ library. But that didn't fix the scalar operator results. I then found a whole collection of scalar x scalar/matrix routines in xpow.cc and then replaced all of those with the new overloaded routine (zzpow ()). Consequently, the casting to int cases I discarded; they seemed out-dated in the sense neither of the pow() functions have an integer exponent version. There is

but the documentation says that most of those end up casting to a double...so it is like casting to an int just to be cast back to a double. Maybe I didn't implement that correctly. Anyway, it was producing errors so I took out the static_cast<int> which simplified some code. See if you agree with what I changed there.

Here are a couple other important cases:

I'm a bit confused by all the combinations in xpow.cc versus the inline matrix implementations in mx-inlines.cc. Do a

and one will see an awful lot of uses of std::pow(). So maybe those need a cursory review to confirm there aren't more x^0 issues, say for sparse and so on.

And then there is this from the oct-inttypes.cc file:

The above is a thorough consideration of all the base=1 and exponent=0 cases. (Is it worth testing for base==1? I mean, in most cases it will not be 1 and the test is just wasting a few cycles because the routine will produce a valid result if base=1.) So, it's not like this issue hasn't been considered. It's just that the std::pow routine has found extensive use throughout. Then there is a separate powf () which may be repetitive and can be replaced by simple use of the overloaded pow () function. Possibly some cruft here?

Lastly, I removed the following comment:

I checked, and I found that two cases in the example given do in fact now match. As for the (-realX)^I.5 comment, i.e., that this scenario should have the real portion forced to 0 (purely imaginary) I'm a bit hesitant to do that. (We should be able to devise a short command to assess if the exponent has a fractional portion of exactly 0.5 if wanted, but the question is if we want that.) My thinking is that the std::pow routine is using an algorithm to compute this value and if we force a certain subset of the values to be something slightly different than what the algorithm generates, it's as if we are creating a discontinuity, algorithm-wise. That is, say, someone computes a series of numbers that converges to a -real value and looks at the results of std::pow(); there might be some subtle discontinue in some residual that catches the observer off guard. To summarize, I'm thinking that if it is algorithmic in nature, probably best to leave it as is, but in the case of 0^0 it is an undefined input that we are simply assigning a more desired value to. The user could certainly introduce logic as I described above to make (-realX)^I.5 purely imaginary.

(file #44562)

Dan Sebald <sebald>
Sat 14 Jul 2018 05:30:27 PM UTC, comment #4:

The answer is yes:
Short answer is complex(0,0) is complex, 0+0j is not

Michael Godfrey <godfrey>
Project Member
Sat 14 Jul 2018 08:55:09 AM UTC, comment #3:

The ^ operator and power() function call end up using the same op_el_pow construction. It's a long series of defines that builds a table of operators for which lookup is done to get the function to call.

I can see in the file mx-inlines.cc that at least in the case of matrices and vectors (may not be the same for scalar operations, don't know), the individual power computations are done with the routines:

These inline routines use the std::pow() function which is a complex function. [There is a comment within that suggests the compiler is left to choose pow(), but the "using" directive might force std::pow(). The other option would be cmath pow(double, double). There is also a FIXME questioning this. Should clear this up.] The reference for std::pow is here (both C98 and C++11):

http://www.cplusplus.com/reference/complex/pow/

and the reference for cmath pow is here:

http://www.cplusplus.com/reference/cmath/pow/

In both cases, the base of 0 and exponent of 0 (either double or complex) results in a returned value which is application specific, whereas the error flags that are thrown are more well-defined. So, it seems to me either the options are to check for the exponent being 0 prior, or check for error flags after calling pow() and then decipher. Probably the former.

This creates a problem, though, in the sense that these inline routines are defined as templates and we can't generally do

The following fails as well for cases where Y and R are not complex, e.g., just "float":

So, there is a not-so-easy task of either expanding out all the different scenarios based on R, X and Y typenames. Or, write the two-operand mx_inline_XYZ routines to have another two inputs, call them ZERO and ONE as follows

where the y is compared against ZERO and if equal then r = ONE. But that means all two-argument operator functions have to change, even though the inlining would completely ignore those last two values if they don't appear within the operator function.

This isn't a fun change.

Dan Sebald <sebald>
Sat 14 Jul 2018 06:12:52 AM UTC, comment #2:

As a curious aside, is there a difference between the following?

Dan Sebald <sebald>
Fri 13 Jul 2018 01:10:10 PM UTC, comment #1:

I confirm that Matlab gives 1 in all the cases.

Marco Caliari <caliari>
Project Member
Fri 13 Jul 2018 07:22:49 AM UTC, original submission:

Normally 0.^0 is 1, but in some case it results in NaN.

Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by sebald (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

     

     

    Follow 2 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-07-16 sebald Attached File- => Added octave-zero_to_zero_power-djs2018jul15.patch, #44562
    2018-07-13 caliari StatusNone => Confirmed

    Back to the top


    Powered by Savane 3.3