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.
|
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.
|
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.
|
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.
|
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.
|
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.
|
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.
|
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.
|
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).
|
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.
|
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.
|
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.
|
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)
|
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)
|
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.
|
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)
|
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?
|
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.
|
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.
|
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).
|
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.
|
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.
|
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.
|
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.
|
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)).
|
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.
|
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.
|
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
|
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
|
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.
|