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

 
 

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

Submitter:  Ricardo Biloti <biloti>
Submitted:  Thu 02 Oct 2014 02:31:48 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Fixed Assigned to:  None
Originator Name:  Ricardo Biloti Open/Closed:  * Closed
Release:  * 3.8.1 Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

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>
Group 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>
Group 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>
Group 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>
Group 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 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>
Group 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 x*x 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 x*x 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>
Group 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


octave:12> asin ([1 2])
ans =

   1.5708 + 0.0000i   1.5708 - 1.3170i

octave:13> asin ([1i 2])
ans =

   0.00000 + 0.88137i   1.57080 + 1.31696i


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.


asin ([2 + i/1e10, 2 - i/1e10])
ans =

   1.5708 + 1.3170i   1.5708 - 1.3170i


Rik <rik5>
Group 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>
Group administrator
Tue 07 Oct 2014 08:58:04 PM UTC, comment #20: 

From Matlab's documentation:


On the branch cut, the values of arcsin coincide with the limit "from below" for real arguments x > 1. The values coincide with the limit "from above" for real x < - 1:


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:


octave:1> asin (complex (2,0))
ans =  1.5708 + 1.3170i
octave:2> asin (complex (2,-0))
ans =  1.5708 - 1.3170i


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

Rik <rik5>
Group 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>
Group 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%.


  return fabs (x) > 1.0 ? asin (Complex (x, x < 0 ? 0.0 : -0.0))
                        : Complex (asin (x));


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>
Group administrator
Tue 07 Oct 2014 05:14:58 AM UTC, comment #17: 

The phase makes a difference because


x*x


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


2-0i * 2-0i = 4 - 0i


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.


 Complex
 rc_asin (double x)
 {
-  return fabs (x) > 1.0 ? asin (Complex (x)) : Complex (asin (x));
+  return fabs (x) > 1.0 ? asin (Complex (x, signbit (x) ? 0.0 : -0.0))
+                        : Complex (asin (x));
 }


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.


rik Patch
---------
2.277
2.261
2.263
2.257

mean = 2.2645

jwe Patch
---------
2.211
2.213
2.217
2.234

mean = 2.2188



(file #32237)

Rik <rik5>
Group 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>
Group 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>
Group 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>
Group 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>
Group 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.


  return fabs (x) > 1.0 ? asin (Complex (x, copysign (x, 0.0)))
                        : Complex (asin (x));


Rik <rik5>
Group 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>
Group 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>
Group 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


template<> class complex<double> {
  complex (double re = 0.0, double im = 0.0);


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


   return fabs (x) > 1.0 ? asin (Complex (x)) : Complex (asin (x));


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.


Complex
rc_asin (double x)
{
  return fabs (x) > 1.0 ? asin (Complex (x, signbit (x) ? -0.0 : +0.0))
                        : Complex (asin (x));
}


Rik <rik5>
Group 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


diff --git a/liboctave/numeric/lo-mappers.cc b/liboctave/numeric/lo-mappers.cc
--- a/liboctave/numeric/lo-mappers.cc
+++ b/liboctave/numeric/lo-mappers.cc
@@ -496,7 +496,17 @@
 Complex
 rc_asin (double x)
 {
-  return fabs (x) > 1.0 ? asin (Complex (x)) : Complex (asin (x));
+  if (fabs (x) > 1.0)
+    {
+      static Complex i (0, 1);
+
+      // Don't forward to asin (Complex (x)) because the value of
+      // Complex (1.0 - x*x) != 1.0 - (Complex (x) * Complex (x))!
+
+      return -i * log (i*x + sqrt (Complex (1.0 - x*x)));
+    }
+  else
+    return Complex (asin (x));
 }

 FloatComplex


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

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

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


sqrt (-3)


but in the C++ code we are computing


sqrt (1.0 - Complex (2, 0) * Complex (2, 0))


which is equivalent to


sqrt (Complex (-3, -0))


Oops.

John W. Eaton <jwe>
Group 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


-i * log (i*x + sqrt (1.0 - x*x))


produces the correct answer


octave:9> -i * log (2i + sqrt (-3))
ans =  1.5708 - 1.3170i
octave:10> -i * log (-2i + sqrt (-3))
ans = -1.5708 + 1.3170i


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

Rik <rik5>
Group 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:


Complex
asin (const Complex& x)
{
  static Complex i (0, 1);

  return -i * log (i*x + sqrt (1.0 - x*x));
}


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


  return -i * log (i*x + signum (real (x)) * sqrt (1.0 - x*x));


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:


On the branch cut, the values of arcsin coincide with the limit "from below" for real arguments x > 1. The values coincide with the limit "from above" for real x < - 1:


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>
Group 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>
Group 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


y = asin (2)



Rik <rik5>
Group 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 Files:
   
   
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
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by biloti (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2014-10-10 rik5 StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2014-10-07 rik5 Attached File- Added asin.patch, #32237
    2014-10-06 jwe Attached File- Added diffs.txt, #32234
    2014-10-06 jwe Attached File- Added diffs.txt, #32233
    2014-10-06 rik5 Summaryasin behaves differently from Matlab for arguments larger than 1 asin() behaves differently from Matlab for arguments larger than 1
    2014-10-06 rik5 StatusNeed Info Confirmed
    2014-10-06 rik5 StatusNone Need Info

    Back to the top

    Powered by Savane 3.13-758e.
    Corresponding source code