bugGNU Octave - Bugs: bug #43349, asin() behaves differently from...

 
 

bug #43349: asin() behaves differently from Matlab for arguments larger than 1

Submitted by:  Ricardo Biloti <biloti>
Submitted on:  Thu 02 Oct 2014 02:31:48 PM UTC  
 
Category: Octave FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Matlab Compatibility
Status: FixedAssigned to: None
Originator Name: Ricardo BilotiOpen/Closed: Closed
Release: 3.8.1Operating 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 Jump to the original submission)

Fri 10 Oct 2014 03:42:49 PM UTC, comment #29:

Thank you all for your attention to this issue. Now it is clear where the problem is and how to solve it, I will implement a workaround to circumvent it until the fix be officially released.

Ricardo Biloti <biloti>
Fri 10 Oct 2014 12:52:12 AM UTC, comment #28:

I checked in a fix to asin and acos based on jwe's suggestion (http://hg.savannah.gnu.org/hgweb/octave/rev/ba7e42dea4b2). I updated the NEWS file to mention this since people might have scripts that were somehow depending on the old behavior.

@Ricardo: The fix was made on the development branch of Octave. In order to access the fix you will need to build from Mercurial sources. Instructions for that are on the Octave web page. Otherwise, this will be a part of the 4.2 release, but I wouldn't expect to see that for maybe a year.

Rik <rik5>
Project Administrator
Thu 09 Oct 2014 02:48:12 PM UTC, comment #27:

This is a minor point, but I do not think
that this is correctly classified as "Matlab Compatibility".
"Correct result" is more accurate and it happens
to be R and Matlab and other compatible.
Correct results and Matlab compatibility are NOT
the same thing.

Michael Godfrey <godfrey>
Project Member
Thu 09 Oct 2014 03:17:07 AM UTC, comment #26:

I've taken your patch as the basis for the changeset. I'm adding some BIST tests to mappers.cc to check the values on either side of the branch cut to make sure we are Matlab compatible and that we don't accidentally undo this cset at a later date.

Rik <rik5>
Project Administrator
Thu 09 Oct 2014 02:42:00 AM UTC, comment #25:

Again? Ugh.

I think my patch https://savannah.gnu.org/support/download.php?file_id=32234 is the right thing to do. If the imaginary part of x is -0, then x*x will have +0 as the imaginary part, and that looks like a real number, so 1-x*x should not generate a negative imaginary part. Likewise if x has +0 as the imaginary part, it looks like a real number, so x*x should be real, so 1-x*x should not have a negative imaginary part. Does that make sense? In any case and FWIW, I think it is compatible with Matlab.

John W. Eaton <jwe>
Project Administrator
Thu 09 Oct 2014 02:40:41 AM UTC, comment #24:

Grr, bitten by the smart markup features of the bug tracker.

Here it is again:

I think my patch https://savannah.gnu.org/support/download.php?file_id=32234 is the right thing to do. If the imaginary part of x is -0, then x*x will have +0 as the imaginary part, and that looks like a real number, so 1-x*x should not generate a negative imaginary part. Likewise if x has +0 as the imaginary part, it looks like a real number, so xx should be real, so 1-x*x should not have a negative imaginary part. Does that make sense? In any case and FWIW, I think it is compatible with Matlab.

John W. Eaton <jwe>
Project Administrator
Thu 09 Oct 2014 02:34:12 AM UTC, comment #23:

I think my patch https://savannah.gnu.org/support/download.php?file_id=32234 is the right thing to do. If the imaginary part of x is -0, then xx will have +0 as the imaginary part, and that looks like a real number, so 1-x should not generate a negative imaginary part. Likewise if x has +0 as the imaginary part, it looks like a real number, so xx should be real, so 1-x should not have a negative imaginary part. Does that make sense? In any case and FWIW, I think it is compatible with Matlab.

John W. Eaton <jwe>
Project Administrator
Wed 08 Oct 2014 11:54:49 PM UTC, comment #22:

Ugh. This problem just keeps growing. If I make the fix in rc_asin which affects real numbers only, then I get a difference between calling the function on a real array versus a complex array. With the fix in place I see

The problem is that when a complex array is created which has real values in it then the imaginary portion is automatically assigned the value +0.

Thus, I'm now thinking that the solution has to be in asin (Complex) and that jwe's second solution is the correct way to go. The phase of a real number can't be adjusted with a call to complex (x, -0), but that is okay. The convention is that the phase has already been determined by the sign of the real part of x. If one really must force a phase then you can move the number off the real axis by adding a small imaginary component.

Rik <rik5>
Project Administrator
Wed 08 Oct 2014 01:08:05 PM UTC, comment #21:

The case you bring up with asin (complex (2, -0)) does matter and we shouldn't break that while fixing some other problem so the solution used in my first patch is probably better than the second.

But I'm also not sure I've clearly understood all the details here, though at one point I thought I did. So I'd say go ahead with whatever solution you think is best. Or I'd be glad to make the changes I posted in my first patch (with the float/double issues fixed).

John W. Eaton <jwe>
Project Administrator
Tue 07 Oct 2014 08:58:04 PM UTC, comment #20:

From Matlab's documentation:

So they are explicit about which side of the branch cut they use and it depends on the value of the real part of the argument x. My patch tests x to see which side they are on and then adjusts the phase accordingly by flipping the sign bit of the imaginary part (which happens to have value zero).

But, I'm still not certain that it is very portable. I could see alternative compilers not necessarily preserving the sign bit of x mult x. Your solution looks more portable. The only thing I do see is that it essentially does automatic narrowing of -0 to +0. Right now I can adjust the phase if I want by calling the asin (Complex) function directly:

But, I'm not convinced that this edge case is all that important.

Rik <rik5>
Project Administrator
Tue 07 Oct 2014 07:16:00 PM UTC, comment #19:

I don't understand why we should forward (x,0) to the complex asin function if x < 0 but (x,-0) if x >= 0. We're trying to define a function of a real-valued argument, so why should we assign a negative imaginary part for that value?

I'm also not sure that checking for a zero imaginary part as I did in my second patch is really the best solution either. It's just slightly less duplicated code.

I'll look at this a bit more, but I think we should strive for the solution that will be easiest to understand later.

John W. Eaton <jwe>
Project Administrator
Tue 07 Oct 2014 05:27:06 AM UTC, comment #18:

If I swap out the signbit call and just use x < 0 in the test then the performance advantage is down to 0.5%.

I'd say it is now just a question of which is more appealing code. I think modifying the phase angle is less invasive, but it might run into differences between compilers and platform differences on Windows.

Rik <rik5>
Project Administrator
Tue 07 Oct 2014 05:14:58 AM UTC, comment #17:

The phase makes a difference because

preserves the signbit on the imaginary portion of the product when it is zero. I tried in C++ and

The sqrt function distinguishes between the two phase angles -0, +0 and so the results come out different for asin (2-0i) and asin (-2+0i). I had the sign reversed in my first attempt but the attached second attempt works fine.

I used the following as a reference for expected values across the branch cut: http://www.mathworks.com/help/symbolic/mupad_ref/arcsin.html#stdlib-invtrigfuncs-ex4

I think either patch will work. I tried timing them and jwe's method is 2% faster, although that is close to being noise. I tried asin on a 3000 x 3000 element matrix with the following results.

(file #32237)

Rik <rik5>
Project Administrator
Mon 06 Oct 2014 10:04:36 PM UTC, comment #16:

There was some double/float confusion in the patch I posted. Try this version instead.

(file #32234)

John W. Eaton <jwe>
Project Administrator
Mon 06 Oct 2014 09:51:25 PM UTC, comment #15:

Well, one thing to do would be to write the BIST tests for extreme input values for the functions in mappers.cc.

The golden results could be taken from online R or wolframalpha so no one would need a license for Matlab.

We could post that we need help with this project in the Short Projects page on the Octave Wiki, but I haven't seen much follow through on the projects there. We could cross-post to Octave Maintainers and the Help list for this project as well.

Test writing is relatively easy, and once in place would give us an idea of whether all the functions in lo-mappers.cc need rework or only a few.

Rik <rik5>
Project Administrator
Mon 06 Oct 2014 09:48:38 PM UTC, comment #14:

I don't see how setting the sign of the imaginary part to be the same as the real part makes sense or solves the problem.

The issue is that when we have real valued X but compute the subexpression 1-X*X using complex arithmetic we get an imaginary part that is -0 when it really should be +0. And that would happen whether or not you explicitly set the imaginary part of Complex (X) to be 0.

Reviewing all the trig functions would be great, but I think the cases that are subject to this bug are just those that have a term like 1.0 - X*X or 1.0 - X in which imag(X) == 0. Thinking about it now, however, I'm wondering whether we should just check imag (X) == 0 and avoid generating a -0 for the imaginary part regardless. Isn't that the right thing to do, even when the type of the argument is Complex but the imaginary part is 0?

Please try the attached patch.

(file #32233)

John W. Eaton <jwe>
Project Administrator
Mon 06 Oct 2014 09:26:47 PM UTC, comment #13:

I think that you got NaN due to the R requirement of
explicitly making the input complex. That is why
I used cos(asin(2+0i). Without the +0i, R says NaNs...

Is there anyone with suitable skills who could do this
sooner than a year from now?

Michael Godfrey <godfrey>
Project Member
Mon 06 Oct 2014 09:04:38 PM UTC, comment #12:

@Michael: Yes, there would need to be a review of most of the mapping functions in lo-mappers.cc to see if this change needs to be done. Maybe we could interest someone in doing that, or maybe it has to wait until the next GSOC round?

I tried to test in R using an online execution engine but the result I get is NaN. I think it's good enough validation that Mathematica and Matlab use the same convention.

Another possible implementation would be to use copysign instead of the tertiary operator `?'. I don't know which would be faster.

Rik <rik5>
Project Administrator
Mon 06 Oct 2014 09:04:32 PM UTC, comment #11:

R says: R version 3.1.1 (2014-07-10)

cos(asin(2+0i))
[1] 0+1.732051i

It requires that the input be complex.
So, looks OK.

Michael Godfrey <godfrey>
Project Member
Mon 06 Oct 2014 08:48:31 PM UTC, comment #10:

Rik,

That looks right. But, this needs to be done
throughout the trig library, right?

Also, it would be good to check R just to know
(I assume that they are "Matlab conmpatible).

Michael Godfrey <godfrey>
Project Member
Mon 06 Oct 2014 08:25:41 PM UTC, comment #9:

I agree that the problem is the conversion of a real to a complex value with zero for the imaginary unit. I think the problem is that we use the default form of the C++ Complex constructor with just a single argument. The Complex constructor is shown below

When it is called in rc_asin with just a single argument then the imaginary unit takes the default which is always +0.0.

However, if I use the two argument form and explicitly set the phase of the imaginary unit then I get the correct answer. The replacement code below works on the examples in this bug report.

Rik <rik5>
Project Administrator
Mon 06 Oct 2014 05:31:31 PM UTC, comment #8:

Some care is needed here. You don't want to just throw away the signbit on the imaginary part of the argument to sqrt if it happens to be zero. So I think the fix belongs in rc_asin functions, not the asin (Complex) functions. Something like

We probably need to review all of these functions that convert real-valued arguments to complex.

John W. Eaton <jwe>
Project Administrator
Mon 06 Oct 2014 05:17:05 PM UTC, comment #7:

Ugh, the difference is that in the interpreter you are writing

but in the C++ code we are computing

which is equivalent to

Oops.

John W. Eaton <jwe>
Project Administrator
Mon 06 Oct 2014 04:53:07 PM UTC, comment #6:

Actually, something else must be going on. Direct calculation within Octave of the expression used in the calculation

produces the correct answer

So after calculation, as it is passed back to the Octave interpreter something is flipping the sign bit.

Rik <rik5>
Project Administrator
Mon 06 Oct 2014 04:50:06 PM UTC, comment #5:

Definitely, this is due to the choice of the branch cut. What I can tell you is how I ran into this issue. I was solving the Zoeppritz equations for elastic reflection coefficient. It turns out the reflection coefficients computed have the imaginary part with the wrong signal. After spending quite a time looking for a bug in my code, I have decided to try on Matlab. To my surprise, the results were correct there.

So, that branch cut has impacts that goes beyond x=sin(asin(x)).

Ricardo Biloti <biloti>
Mon 06 Oct 2014 04:27:50 PM UTC, comment #4:

The code in question is in liboctave/numeric/lo-mappers.cc:

It would be possible to put a test on the sign of the real part of x to get the same result as Matlab. Conceptually, this is

But I'm not sure it is worth it. sin (asin (x)) returns x in either case which seems to be the more important thing.

According to Matlab's documentation:

Not everyone gets a vote, but I see that Mathematica also returns the same branch cut values as Matlab.

If this gets changed, then probably all of the functions in lo-mappers.cc that involve branch cuts should be reviewed.

Rik <rik5>
Project Administrator
Mon 06 Oct 2014 12:57:05 PM UTC, comment #3:

In Octave or Matlab:

plot(imag(asin(-2:0.1:2)))

shows the problem graphically.

Michael Godfrey <godfrey>
Project Member
Mon 06 Oct 2014 11:58:18 AM UTC, comment #2:

Below I quote some examples comparing Octave (3.8.1) with Matlab 8.0.0.783 (R2012b). I also have access to Matlab 7.8.0.347 (R2009a), which behaves just like the newer Matlab.

Octave
======
octave:1> version
ans = 3.8.1
octave:2> asin([-2 -1 1 2]')
ans =

-1.5708 + 1.3170i
-1.5708 + 0.0000i
1.5708 + 0.0000i
1.5708 + 1.3170i

octave:3> sin(asin([-2 -1 1 2]'))
ans =

-2.0000 + 0.0000i
-1.0000 + 0.0000i
1.0000 + 0.0000i
2.0000 + 0.0000i

octave:4> cos(asin([-2 -1 1 2]'))
ans =

0.00000 + 1.73205i
0.00000 + 0.00000i
0.00000 - 0.00000i
0.00000 - 1.73205i

Matlab
======

>> version


ans =

8.0.0.783 (R2012b)

>> asin([-2 -1 1 2]')


ans =

-1.5708 + 1.3170i
-1.5708
1.5708
1.5708 - 1.3170i

>> sin(asin([-2 -1 1 2]'))


ans =

-2.0000 + 0.0000i
-1.0000
1.0000
2.0000 - 0.0000i

>> cos(asin([-2 -1 1 2]'))


ans =

0.0000 + 1.7321i
0.0000
0.0000
0.0000 + 1.7321i

Ricardo Biloti <biloti>
Mon 06 Oct 2014 04:58:42 AM UTC, comment #1:

I think I see what you mean. To make it concrete, could you report the version of Matlab you are using and also what a single specific example returns, such as

Rik <rik5>
Project Administrator
Thu 02 Oct 2014 02:31:48 PM UTC, original submission:

Let x be a real value greater than 1 (x > 1). Then Matlab reports that

asin (x) = pi/2 - a*i, for some a > 0

However Octave produces the complex conjugate

asin (x) = pi/2 + a*i

This occurs only for real x > 1. Both Matlab and Octave agree for real x <=1.

Ricardo Biloti <biloti>

 

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

Attach File(s):
   
   
Comment:
   

Attached Files
file #32237:  asin.patch added by rik5 (704B - text/x-diff)
file #32234:  diffs.txt added by jwe (2KiB - text/plain)
file #32233:  diffs.txt added by jwe (2KiB - text/plain)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by jwe (Posted a comment)
  • -unavailable- added by godfrey (Posted a comment)
  • -unavailable- added by rik5 (Posted a comment)
  • -unavailable- added by biloti (Submitted the item)
  •  

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

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 8 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Fri 10 Oct 2014 12:52:12 AM UTCrik5StatusConfirmed=>Fixed
      Open/ClosedOpen=>Closed
    Tue 07 Oct 2014 05:14:58 AM UTCrik5Attached File-=>Added asin.patch, #32237
    Mon 06 Oct 2014 10:04:36 PM UTCjweAttached File-=>Added diffs.txt, #32234
    Mon 06 Oct 2014 09:48:38 PM UTCjweAttached File-=>Added diffs.txt, #32233
    Mon 06 Oct 2014 08:26:05 PM UTCrik5Summaryasin behaves differently from Matlab for arguments larger than 1=>asin() behaves differently from Matlab for arguments larger than 1
    Mon 06 Oct 2014 04:27:50 PM UTCrik5StatusNeed Info=>Confirmed
    Mon 06 Oct 2014 04:58:42 AM UTCrik5StatusNone=>Need Info

    Back to the top


    Powered by Savane 3.1-cleanup1