bugGNU Octave - Bugs: bug #45507, acos returns different results on...

 
 

bug #45507: acos returns different results on big input values on gcc 4.8

Submitter:  Hartmut <hardy>
Submitted:  Thu 09 Jul 2015 06:53:00 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Regression
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Tue 19 Sep 2017 07:25:54 PM UTC, comment #58: 

I think rebasing mxe-octave on msys2 might bring in a new set of system libraries. I have no idea whether this would include an up-to-date libm, or even whether this is considered fixable or if the error is fundamental to the Windows math libraries.

I don't think we want to return this test to an "xtest", because we don't want it to just fail and say "expected failure" on GNU/Linux systems where this must work correctly.

We do not have a way to flag a test "if this fails it is a regression, but only on this operating system, on this other one it's ok to be completely wrong".

Is anyone building Octave natively on Windows using mingw / msys or msys2? I think it would be valuable if some of our Windows users following the default branch would set up a real msys or msys2 development environment on Windows.

Mike Miller <mtmiller>
Group Member
Mon 18 Sep 2017 03:58:41 PM UTC, comment #57: 

Found an old build from July 2016 (hg id 289409b2992d) which failed on Windows for the same tests. So it looks like inverse trigonometric function did not work for large values in Windows since a very long time (ever?).
The comments above these tests suggest to upgrade libm. Is that something that we can do in mxe-octave?

Markus Mützel <mmuetzel>
Group administrator
Sun 17 Sep 2017 08:36:05 PM UTC, comment #56: 

These tests now appear as regressions on Windows 10 64bit:

D:\Repositories\Octave\test\octave-2017-09-17-18-08_clean_gcc7_2_qt5\octave-2017-09-17-18-08\share\octave\4.3.0+\etc\tests\libinterp\corefcn\mappers.cc-tst
***** test <*45507>
 x = [1, -1, i, -i] .* 1e150;
 v = [0, pi, pi/2, pi/2];
 assert (real (acos (x)), v);
!!!!! regression: http://octave.org/testfailure/?45507
ASSERT errors for:  assert (real (acos (x)),v)

  Location  |  Observed  |  Expected  |  Reason
    (1)         1.5708         0         Abs err 1.5708 exceeds tol 0
    (2)         1.5708       3.1416      Abs err 1.5708 exceeds tol 0
shared variables   scalar structure containing the fields:

    rt2 =  1.4142
    rt3 =  1.7321
[snip]
***** test <*45507>
 x = [1, -1, i, -i] .* 1e150;
 v = [0, pi, pi/2, -pi/2];
 assert (imag (acosh (x)), v);
!!!!! regression: http://octave.org/testfailure/?45507
ASSERT errors for:  assert (imag (acosh (x)),v)

  Location  |  Observed  |  Expected  |  Reason
    (4)           -0        -1.5708      Abs err 1.5708 exceeds tol 0
shared variables   scalar structure containing the fields:

    rt2 =  1.4142
    rt3 =  1.7321
***** test <*45507>
 x = [1, -1, i, -i] .* 1e150;
 v = [pi/2, -pi/2, 0, -0];
 assert (real (asin (x)), v);
!!!!! regression: http://octave.org/testfailure/?45507
ASSERT errors for:  assert (real (asin (x)),v)

  Location  |  Observed  |  Expected  |  Reason
    (1)           0          1.5708      Abs err 1.5708 exceeds tol 0
    (2)           0         -1.5708      Abs err 1.5708 exceeds tol 0
shared variables   scalar structure containing the fields:

    rt2 =  1.4142
    rt3 =  1.7321
***** test <*45507>
 x = [1, -1, i, -i] .* 1e150;
 v = [0, 0, pi/2, -pi/2];
 assert (imag (asinh (x)), v);
!!!!! regression: http://octave.org/testfailure/?45507
ASSERT errors for:  assert (imag (asinh (x)),v)

  Location  |  Observed  |  Expected  |  Reason
    (4)           0         -1.5708      Abs err 1.5708 exceeds tol 0
shared variables   scalar structure containing the fields:

    rt2 =  1.4142
    rt3 =  1.7321


This might be related to the other failing tests in bug #49091. However, those are accuracy related and these tests are off by pi/2. So it might be a different bug in mingw-w64.

I looked through a few older test runs on Windows (during the last 6 months) and these tests seem to have always failed. Should these be marked as xtest again?

Markus Mützel <mmuetzel>
Group administrator
Sat 02 Sep 2017 12:04:27 PM UTC, comment #55: 

Changed to regular BIST tests and marked as regressions here (http://hg.savannah.gnu.org/hgweb/octave/rev/d511f7e15652).  Closing report.

Rik <rik5>
Group administrator
Sat 02 Sep 2017 12:24:37 AM UTC, comment #54: 

I think it would be best to rewrite the four "%!xtest <45507>" as "%!test <*45507>" and close this bug as fixed again. If these tests fail they may show up as regressions in someone's test log, but they will be directed to this bug report, and there is plenty explaining here why these functions may be less accurate on systems with a buggy glibc.

Mike Miller <mtmiller>
Group Member
Sat 02 Sep 2017 12:07:48 AM UTC, comment #53: 

I believe that the C++11 conversion has been done.
Are there reasonable tests in place? It appears so, but
just asking...

In any case good to close this since much of it is no
longer relevant.

Michael Godfrey <godfrey>
Group Member
Fri 01 Sep 2017 11:26:08 PM UTC, comment #52: 

Can we close this bug report?  It seems like this is upstream in glibc and only on old systems;  This is not necessarily an Octave problem.

Rik <rik5>
Group administrator
Mon 04 Jul 2016 02:09:39 AM UTC, comment #51: 

Great.  "test" sounds fine then.

LMDE used to go with Debian Testing, but now goes with Debian Stable.

I think that in the process, my apt-sources have got messed up and so I'm frozen with a system a few years old.  I plan to upgrade once I get a few contiguous hours...

Lachlan Andrew <lachlan>
Sun 03 Jul 2016 11:37:03 PM UTC, comment #50: 

Yes I'm fine with setting this bug aside for now, either with or without the xtests, and moving on to more critical things.

Mike Miller <mtmiller>
Group Member
Sun 03 Jul 2016 11:36:30 PM UTC, comment #49: 

Ok, thanks, that clears up a bit.

I think it may have been this bug (just by searching the commit log for suspicious entries)

https://sourceware.org/bugzilla/show_bug.cgi?id=14996

that was fixed by this commit

https://sourceware.org/git/?p=glibc.git;a=commitdiff;h=a9708fed

in glibc version 2.18.

If this is right, affected systems include Debian 7 ("oldstable" at this point), Ubuntu 12.04 (ancient LTS at this point), CentOS/RHEL 6 (old stable RHEL at this point), CentOS/RHEL 7 surprisingly (current stable RHEL), and your version of LMDE (which I don't know how to check or search for package lists).

If LMDE is supposed to be based on Debian, I'm surprised that it has older packages than Debian 8. And I thought it was supposed to roll with Debian testing, but maybe that has changed.

Mike Miller <mtmiller>
Group Member
Sun 03 Jul 2016 11:35:04 PM UTC, comment #48: 

I forgot to add:

How about we leave this aside ("in progress" and with a failing test if that is what JWE chooses) until all of the crash bugs are fixed?

Lachlan Andrew <lachlan>
Sun 03 Jul 2016 10:35:29 PM UTC, comment #47: 


cc     c99-acos-test.c   -o c99-acos-test -lm
z = 1.5708 + -infi
error: fails assertion real(z) == 0



g++ -std=gnu++11    cxx11-acos-test.cc   -o cxx11-acos-test
./cxx11-acos-test
z = 1.5708 + -infi
error: fails assertion real(z) == 0


OS: Linux Mint Debian Edition
gcc version 4.8.2 (Debian 4.8.2-1)
GNU libc version: 2.17
x86_64

Lachlan Andrew <lachlan>
Sun 03 Jul 2016 07:44:18 PM UTC, comment #46: 

Attached for real this time…

(file #37692, file #37693)

Mike Miller <mtmiller>
Group Member
Sun 03 Jul 2016 07:41:35 PM UTC, comment #45: 

I would start by going back to basics and making sure what we're seeing with Octave is reflective of the system libraries, versus some other effect going on.

The attached programs use the standard library routines to compute the assertion that is failing in Octave, specifically the single assertion "assert(real(acos(1e150)), 0)" that seems to be failing for you and others. One is a C99 program and one is a C++11 program, both should produce identical results.

What confuses me even more is I now have a build of Octave with gcc 4.8.4 on a Ubuntu 14.04 VM, and I don't see a failure in Octave, and both programs succeed for me. So now I'm not sure on which systems we see success or failure. Can we get some details from those of you who see failing tests: OS/distribution version, gcc version, libc version, and architecture?

Mike Miller <mtmiller>
Group Member
Sun 03 Jul 2016 07:13:10 AM UTC, comment #44: 

John, I've been convinced that we shouldn't get rid of the tests entirely.  I like Rik's suggestion of "xtest" for the release, and "test" again afterwards.  I've posted my reasoning to the mailing list.

Mike, that's good sleuthing.  Are you suggesting taking the code from a recent version of glibc and using that if the system library suffers from roundoff?  That sounds like a good option for after the 4.2.0 rush.

In light of JWE's comments, I'm re-opening this as "In progress".

Lachlan Andrew <lachlan>
Sat 02 Jul 2016 04:14:49 PM UTC, comment #43: 

FTR I think the fixes to these functions occurred in libm somewhere between glibc 2.19 and 2.22 (because we have Ubuntu 14.04 users saying it fails, but it works for me on Debian stretch). I don't think the version of GCC matters for this, since the C++ functions just forward the arguments to the appropriate libc functions.

So some kind of configure test could be done for these calls, but we'd need someone to contribute some kind of wrapper, possibly pulled from glibc, that would stand in for the broken function(s).

Mike Miller <mtmiller>
Group Member
Sat 02 Jul 2016 03:52:36 PM UTC, comment #42: 

Lachlan:  I don't think anyone likes to see tests that fail, but removing the tests is definitely not what I'd like to see...

I'm not even sure I like to see xtests because then people just get used to the idea that the tests fail, that they are expected to fail, and then don't even look to see why or whether they are problems that should really be fixed.

If there are ways that we could fix these problems on systems where the functions return incorrect results, then we could use wrappers to do that.  That's really not much different from what we were doing with gnulib since we needed to tag all replaced functions with gnulib:: anyway.

John W. Eaton <jwe>
Group administrator
Thu 30 Jun 2016 09:14:52 PM UTC, comment #41: 

I made the 4 failing tests into xtests:

http://hg.savannah.gnu.org/hgweb/octave/rev/18af3ea1d387

It appears that some systems may have some tolerance errors also, I did not address those in this change.

If someone wants to, some eps tolerances can be added to the failing tests, but we have released before with failing tests on older systems with older versions of the system math libraries.

Mike Miller <mtmiller>
Group Member
Thu 23 Jun 2016 03:13:09 PM UTC, comment #40: 

Another possibility is to use "%!xtest" instead of "%!test" for one more version of Octave.

Rik <rik5>
Group administrator
Thu 23 Jun 2016 02:46:45 PM UTC, comment #39: 

Oops, I meant the "tests" not be dropped.

Mike Miller <mtmiller>
Group Member
Thu 23 Jun 2016 02:46:23 PM UTC, comment #38: 

I would prefer that the functions not be dropped, since they are testing exactly the cases that this change was meant to address.

I don't think we have a "policy" for tests that are known to fail on some systems because of unresolvable problems. They are just known to fail, users see them every time they build and learn to ignore them. Check the mailing list archives for recent 4.0 releases and see the test suite results that people post.

I used to regularly build Octave on RHEL 5 and always get failing tests due to tolerance problems on functions like log, pow, and sqrt, and I learned to ignore them.

Mike Miller <mtmiller>
Group Member
Thu 23 Jun 2016 02:22:18 PM UTC, comment #37: 

I agree they're corner cases.

If there is no change to the code, I'd prefer the failing tests to be dropped.  I don't like "make check" failing.  I also don't like forcing upgrades on users.  (There's already too much e-waste...)

Is there a way to skip the tests based on the gcc version?  That would also solve the problem with xpow and realpow failing for 4.8.

Lachlan Andrew <lachlan>
Thu 23 Jun 2016 01:42:58 PM UTC, comment #36: 

These are fairly extreme corner cases.  My feeling is that we shouldn't spend more valuable programmer time on this issue, when there are more pressing ones.

Eventually, probably within the next 12 months, 4.2 will be released making C++ a requirement and most people are going to meet that with a modern version of gcc.

Rik <rik5>
Group administrator
Thu 23 Jun 2016 11:05:43 AM UTC, comment #35: 

Yes, the #else code also fails.

I was thinking of having an #else clause along the lines of file #37293 (i.e., Marco's code translated to C++).

This is marked as a regression.  I'd suggest we use the code from 3.8.2, but that seems to be as prone to overflow as the current #else code.  (Compiling 3.8.2 fails with my version of bison, so I can't check.)

It's not clear if this remaining aspect of the bug is worth fixing.  It's up to you.

Lachlan Andrew <lachlan>
Thu 23 Jun 2016 03:00:59 AM UTC, comment #34: 

If I understand the history of this bug report, the original bug was the the existing acos, etc, functions overflow when given large magnitude arguments. The original set of tests were all passing, but those tests are only on arguments with small magnitudes. Then you proposed a patch to fix this, and Marco suggested using standard Fortran library routines instead, and I showed that the C++ equivalents work (given a new enough GCC).

Based on that, I would think that the current #else condition will also fail to give correct answers for the new set of tests that I added when building with GCC 4.8. Can you test that by fudging the config.h?

So do we need an entirely new code branch anyway for the #else case if we want to support systems with the C++ functions but bad results from those functions?

Mike Miller <mtmiller>
Group Member
Thu 23 Jun 2016 12:23:28 AM UTC, comment #33: 

Mike, sorry for the slow reply.

Yes, I still get bad results with the patch applied.

Configure says "yes" for std::acos in <complex> (and all the others).

I'm reopening this, because now there are three new failing test on my system (in addition to the two failures from xpow and realpow).  GCC 4.8 isn't really old; as Julien Bect pointed out, even gcc 4.7 is used in ubuntu 12.04, which is not yet end-of-life.

One option would be to have an #if defined () to use std::acos on systems where it works, and the slow explicit code where it doesn't, with a FIXME comment to say it should be removed once 4.8 is out of circulation.


 rval = 0;
 ival = 1.31695789692481635;
 obs = acos ([2, 2-i*eps, 2+i*eps]);
 exp = [rval + ival*i, rval + ival*i, rval - ival*i];
 assert (obs, exp, 2*eps);
 rval = pi;
 obs = acos ([-2, -2-i*eps, -2+i*eps]);
 exp = [rval - ival*i, rval + ival*i, rval - ival*i];
 assert (obs, exp, 2*eps);
 assert (acos ([2 0]),  [ival*i, pi/2], 2*eps);
 assert (acos ([2 0i]), [ival*i, pi/2], 2*eps);
 x = [1, -1, i, -i] .* 1e150;
 v = [0, pi, pi/2, pi/2];
 assert (real (acos (x)), v);
!!!!! test failed
ASSERT errors for:  assert (real (acos (x)),v)

  Location  |  Observed  |  Expected  |  Reason
    (1)         1.5708         0         Abs err 1.5708 exceeds tol 0
    (2)         1.5708       3.1416      Abs err 1.5708 exceeds tol 0
shared variables
  scalar structure containing the fields:

    rt2 =  1.4142
    rt3 =  1.7321
***** test
 rval = pi/2;
 ival = 1.31695789692481635;
 obs = asin ([2, 2-i*eps, 2+i*eps]);
 exp = [rval - ival*i, rval - ival*i, rval + ival*i];
 assert (obs, exp, 2*eps);
 obs = asin ([-2, -2-i*eps, -2+i*eps]);
 exp = [-rval + ival*i, -rval - ival*i, -rval + ival*i];
 assert (obs, exp, 2*eps);
 assert (asin ([2 0]),  [rval - ival*i, 0], 2*eps);
 assert (asin ([2 0i]), [rval - ival*i, 0], 2*eps);
 x = [1, -1, i, -i] .* 1e150;
 v = [pi/2, -pi/2, 0, -0];
 assert (real (asin (x)), v);
!!!!! test failed
ASSERT errors for:  assert (real (asin (x)),v)

  Location  |  Observed  |  Expected  |  Reason
    (1)           0          1.5708      Abs err 1.5708 exceeds tol 0
    (2)           -0        -1.5708      Abs err 1.5708 exceeds tol 0
shared variables
  scalar structure containing the fields:

    rt2 =  1.4142
    rt3 =  1.7321
***** test
 v = single ([0, pi/2*i, 0, -pi/2*i]);
 x = single ([0, i, 0, -i]);
 assert (asinh (x), v,  sqrt (eps ("single")));
 x = [1, -1, i, -i] .* 1e150;
 v = [0, 0, pi/2, -pi/2];
 assert (imag (asinh (x)), v);
!!!!! test failed
ASSERT errors for:  assert (imag (asinh (x)),v)

  Location  |  Observed  |  Expected  |  Reason
    (4)           -0        -1.5708      Abs err 1.5708 exceeds tol 0
shared variables
  scalar structure containing the fields:

    rt2 =  1.4142
    rt3 =  1.7321


Lachlan Andrew <lachlan>
Fri 17 Jun 2016 09:55:34 PM UTC, comment #32: 

Fixed on the default branch now. In the absence of any test suggestions, I added some minimal unit testing of large magnitude arguments.

http://hg.savannah.gnu.org/hgweb/octave/rev/e9765b62d4e8

Mike Miller <mtmiller>
Group Member
Wed 15 Jun 2016 01:30:34 AM UTC, comment #31: 

Lachlan: one more update to the patch attached. So you're saying with this patch applied, you still get bad results with gcc 4.8? Does configure say yes or no to "checking for std::acos in <complex>..."?

Marco: thanks for the fix, now all values look correct.

(file #37480)

Mike Miller <mtmiller>
Group Member
Tue 14 Jun 2016 11:38:15 PM UTC, comment #30: 

So is the plan to leave this bug unfixed for gcc 4.8.2?

We could use compiler-dependent code to use std::acos etc. for more recent versions and the explicit code for gcc 4.8.2, or just say that 4.8.2 is limited to small input values.

Lachlan Andrew <lachlan>
Tue 14 Jun 2016 07:01:06 PM UTC, comment #29: 

asinh(1e8j) was my fault in textx3. Please use testx4, in which I simplified acosh cases, too.

(file #37476)

Marco Caliari <caliari>
Group Member
Tue 14 Jun 2016 03:38:10 PM UTC, comment #28: 

Thanks, I combine that with my empirical testing to mean that the result of std::acos should be conjugated only when z is purely real and positive. Updated patch attached, also covering bug #44310.

Almost all test cases in Marco's testx3 script now match. Only mismatching one is


>> asinh (1e8j)
ans =  19.1138 +  1.5708i


So is the postprocessing here to invert the sign of the real part when real(z) == 0 and imag(z) > 1?

I included std::atan as well for symmetry, even though it may not have had any domain problems as written.

Anyone care to contribute any test cases to this, e.g. based on testx3?

(file #37472)

Mike Miller <mtmiller>
Group Member
Tue 14 Jun 2016 06:27:30 AM UTC, comment #27: 

Hi Mike,

from the definition here

http://it.mathworks.com/help/matlab/ref/acos.html

isn't it just negative imaginary part if abs(z)>1, z<0 and positive imaginary part if abs(z)>1, z>0?

Marco Caliari <caliari>
Group Member
Mon 13 Jun 2016 11:36:32 PM UTC, comment #26: 

Hmm, I get this ("ans" lines removed):


acos
    0.000 + 19.114i    0.000 + 19.114i
    3.142 - 19.114i    3.142 + 19.114i
    1.571 - 19.114i    1.571 - 19.114i
    1.571 + 19.114i    1.571 + 19.114i
    0.785 - 19.460i    0.785 - 19.460i
    0.785 + 19.460i    0.785 + 19.460i
    2.356 - 19.460i    2.356 - 19.460i
    2.356 + 19.460i    2.356 + 19.460i


This looks correct except for one case, so as before, I just need some hint from someone who understands what should be done about the imaginary part to make the result Matlab compatible.


>> acos (-1e8)
ans =   3.1416 + 19.1138i


According to Marco's test script this should be conjugated, or should not have been conjugated.

Mike Miller <mtmiller>
Group Member
Sat 28 May 2016 07:35:57 AM UTC, comment #25: 

Mike, with gcc 4.8.2 and your patch (updated to reflect the addition of octave namespace -- see attached), Marco's test from comment #3 gives


acos
ans =
    0.000 + 19.114i    1.571 +    Infi
ans =
    3.142 - 19.114i    1.571 +    Infi
ans =
    1.571 - 19.114i    1.571 -    Infi
ans =
    1.571 + 19.114i    1.571 + 19.114i
ans =
    0.785 - 19.460i    1.571 -    Infi
ans =
    0.785 + 19.460i    0.785 + 19.460i
ans =
    2.356 - 19.460i    1.571 -    Infi
ans =
    2.356 + 19.460i    2.356 + 19.460i


(file #37293)

Lachlan Andrew <lachlan>
Thu 26 May 2016 03:40:24 AM UTC, comment #24: 

Ping. My patch is ready, uses the std::acos library function that operates on std::complex variables, just needs some touching up from someone who understands why the Matlab acos imaginary part differs from the standard library.

Mike Miller <mtmiller>
Group Member
Fri 04 Mar 2016 06:51:59 PM UTC, comment #23: 

My patch still needs work, fails a test in make check:


>> acos (-2)
ans =  3.14159265358979 + 1.31695789692482i


I only guessed at how to correct the phase by looking at the current code, I am no expert at complex analysis. If someone can tell me how to fix the sign of the output of std::acos so it corresponds to Matlab and other numerical software I can incorporate that.

Mike Miller <mtmiller>
Group Member
Fri 04 Mar 2016 01:13:43 AM UTC, comment #22: 

Mike, don't worry about throwing out my attempts.  I had assumed from Comment #13 that they would be discarded.  Your approach looks much nicer!

Lachlan Andrew <lachlan>
Fri 04 Mar 2016 12:32:01 AM UTC, comment #21: 

Is something like the attached what we have in mind for this bug? I'm coming in late so I'm not following all of the discussion. And I don't want to throw out Lachlan's hard work on this, but it seemed like the discussion turned to "can't we do this a simpler way?" This patch uses the std::acos routine provided by C++, and corrects the phase if I'm understanding that part correctly.

(file #36541)

Mike Miller <mtmiller>
Group Member
Wed 10 Feb 2016 03:31:05 AM UTC, comment #20: 

John,

This seems reasonable for now. It might be useful
to include a // FIXME comment about the current need
to support systems that lack required math libraries.

Michael Godfrey <godfrey>
Group Member
Tue 09 Feb 2016 06:42:01 AM UTC, comment #19: 

Mike, I agree that we should avoid mixing C and C++ complex is probably not desirable, but in this case it might be better than any other alternative.

We could write our function something like this:


std::complex<double>
xacos (const std::complex<double>& z)
{
#if defined (HAVE_COMPLEX_STD_ACOS)
  return std::acos (z);
#elif (HAVE_CACOS)
  _Complex double r_c99 = cacos (real (z) + imag (z) * I);

  return std::complex<double> (creal (r_c99), cimag (r_c99));
#else
   // Our best attempt?
   // We could just copy the code from the GNU C library...
#endif
}


John W. Eaton <jwe>
Group administrator
Thu 04 Feb 2016 08:46:22 PM UTC, comment #18: 

I was not aware that libm was dependent on C++11. It is
used at many places in Octave already (possibly conditionally
on various things).

Separately, is it known what systems do not yet provide C++11?

Michael Godfrey <godfrey>
Group Member
Thu 04 Feb 2016 08:28:29 PM UTC, comment #17: 

Octave already prefers std::complex to C99's complex type. Using both would be awkward and unpreferable.

C++11 does provide std::acos(std::complex), but we don't require C++11 at this point and can't guarantee that all platforms Octave currently supports have C++11 available.

Mike Miller <mtmiller>
Group Member
Thu 04 Feb 2016 08:09:04 PM UTC, comment #16: 

Why not use libm?  include complex.h
includes the cacos or variants for double of float
if needed which take complex arg.
See:

https://sourceware.org/newlib/libm.html#cacos

Michael Godfrey <godfrey>
Group Member
Thu 04 Feb 2016 07:09:05 PM UTC, comment #15: 

@Marco: It might be possible to use Fortran libraries, but jwe would need to sign off on that.  In general, we try to link against more modern code like libc, but libc doesn't have an acos function which accepts complex inputs or any inputs outside the domain [-1,1].

Rik <rik5>
Group administrator
Thu 04 Feb 2016 09:40:44 AM UTC, comment #14: 

Hi,

I tried


write(6,*) acos((1.0d8,0.0d0))


in fortran and got (0.0000000000000000,-19.113827924512311). A part from the sign of the imaginary part (which can be easly changed), no overflow due to catastrophic cancellation. Why is it not possible to use the same standard math library used by fortran (gfortran-4.8 in my case) and fix the sign later?

Marco Caliari <caliari>
Group Member
Wed 03 Feb 2016 09:58:17 PM UTC, comment #13: 

The acos routine for a single numeric value is called numel() times when acos is calculated over an array.  It seems that the function has grown significantly more complex (no pun intended).  This may only be appearances because the single math expression is now calculated in distinct sections with flow control between the sub-calculations.  Regardless, one thing to do is benchmark Octave performance before and after the patch.  As a quick test, I used


x = 1e20*rand (1e3,1e3,20) - 5e19;
save -binary bm.var x
tic; y = acos (x); toc


I saved the test values so I could use the exact same numbers before and after the patch was applied.  I ran the tic/toc test 4 times and discarded the most extreme value to leave 3 samples.


mean execution time pre-patch  = 1.4642
mean execution time post-patch = 2.2985
percentage change = +56%


This is obviously going through a slow path in the code because if I use


x = rand (1e3,1e3,15);


and re-test then the results are equivalent.

Is there a way to improve the performance of the slow path?  Or does it just not matter all that much because extremely large inputs to acos may occur, but are rare events?

Rik <rik5>
Group administrator
Wed 03 Feb 2016 06:06:15 PM UTC, comment #12: 

Complex is just a typedef for std::complex<double>, so I think it might be better to define the templates in terms of the real type only:


template <typename T>
std::complex<T>
acos (const std::complex<T>& x)
{
  ...
}


Then could we eliminate the wrapper function?

Please also use "typename" instead of "class" in the template declarations.


John W. Eaton <jwe>
Group administrator
Sat 30 Jan 2016 10:18:58 AM UTC, comment #11: 

Ah, float numbers. So, sqrt(1-x*x)==sqrt(-x*x) if real(x) >= 2^27 or imag(x) >= 2^27. The same is true with 2^13 if x is single type. So, 2 to the power ceil(log2(1/eps)/2+1).

Marco

Marco Caliari <caliari>
Group Member
Sat 30 Jan 2016 09:54:31 AM UTC, comment #10: 

Marco, you're right that the limit is imprecise.  It isn't intended to be the precise value for which the "1+" becomes irrelevant.  It was just intended to be something in the wide range such that it is harmless to skip the sqrts for a double and it will avoid overflow for a float.

I'll use 2^27 in the next version of the patch -- but I'll wait to see if the eventual commiter has any comments before rolling the next version.

Lachlan Andrew <lachlan>
Sat 30 Jan 2016 08:02:52 AM UTC, comment #9: 

@ Andrew

Thanks Andrew for converting my script into a working patch! Nice your trick to avoid sqrt * sqrt. If I'm right, you check whether sqrt(one - x*x) == sqrt(-x*x) for which the square root is trivial. But your limit 1e18 is a bit imprecise. Isn't it enough that fabs (re) > 2^27 or fabs (im) > 2^27?

Marco

Marco Caliari <caliari>
Group Member
Fri 29 Jan 2016 08:25:07 PM UTC, comment #8: 

Thanks for this second patch, Andrew.

I've compiled the patch with the current "default" branch. The results from my original acos calls are now fully Matlab compatible:

acos(1e7) = 0 + 16.81124i
acos(1e8) = 0 + 19.11383i

I can also confirm that "make check" doesn throw any errors, now.



Hartmut <hardy>
Fri 29 Jan 2016 10:11:42 AM UTC, comment #7: 

Thanks for checking, Hartmut.

The problem with the previous patch was the use of 'abs' instead of 'fabs'.

The attached patch passes Marco's tests.  Thanks Marco! That script was really helpful, both as a benchmark and showing how to avoid the cancellations.

It also passes  make check  (although it slightly widens the tolerance on one test).

I've tried to avoid the double sqrt; the patch could be made slightly shorter and presumably slightly slower by using the double square root.  That is a matter of taste for whoever applies the patch.

(file #36198)

Lachlan Andrew <lachlan>
Tue 26 Jan 2016 12:13:07 PM UTC, comment #6: 

I also observe a failing test with his patch (from comment #4):

"test mappers.cc-tst" returns:


***** test
 rval = 0;
 ival = 1.31695789692481635;
 obs = acos ([2, 2-i*eps, 2+i*eps]);
 exp = [rval + ival*i, rval + ival*i, rval - ival*i];
 assert (obs, exp, 2*eps);
 rval = pi;
 obs = acos ([-2, -2-i*eps, -2+i*eps]);
 exp = [rval - ival*i, rval + ival*i, rval - ival*i];
 assert (obs, exp, 2*eps);
 assert (acos ([2 0]),  [ival*i, pi/2], 2*eps);
 assert (acos ([2 0i]), [ival*i, pi/2], 2*eps);
!!!!! test failed
ASSERT errors for:  assert (obs,exp,2 * eps)

  Location  |  Observed  |  Expected  |  Reason
    (1)          0+1i       0+1.317i     Abs err 0.31696 exceeds tol 4.4409e-16
shared variables
  scalar structure containing the fields:

    rt2 =  1.4142
    rt3 =  1.7321


Hartmut <hardy>
Tue 26 Jan 2016 11:49:54 AM UTC, comment #5: 

I have tested the current default branch with this patch (from comment #4). Here are my results for the acos function:

acos(1e7) = 0 + 16i
acos(1e8) = 0 + 19i

So the results still have the same sign in the imaginary part as Matlab (and as Octave 4.0.0 release) = good.

The results don't have the jump to NaN values somewhere between 1e7 and 1e8 any more, which is similar to what Matlab does (and what Octave 3.8.2 release did) = good.

But the results are still slightly different to Matlab: The imaginary part seems to always be INTEGER valued, and not floating point valued. Is this intended?


Hartmut <hardy>
Tue 26 Jan 2016 08:35:20 AM UTC, comment #4: 

The attached patch fixes acos.

I didn't notice Marco's tests; I'll address them all next.

(file #36167)

Lachlan Andrew <lachlan>
Thu 16 Jul 2015 10:04:59 AM UTC, comment #3: 

Hi, I'm on 3.8.2, and I see several discrepancies or NaN (try, e.g., asin(1e8)). Some of them are related to the evaluation of quantities like z+sqrt(z^2-1) which, for large negative numbers z, are z+|z| and therefore 0. Moreover, sqrt(z^2-1) may overflow for z>sqrt(realmax). So, I wrote a test file for acos, asin, acosh and asinh. If I run it on Matlab 2014, the results are always correct. The test file should help in the correct implementation.

Marco

(file #34451)

Marco Caliari <caliari>
Group Member
Fri 10 Jul 2015 07:06:11 AM UTC, comment #2: 

Sorry, I didn't pay enough attention to the sign of the imaginary part. You were right Rik. Here are the double checked results:

Matlab (R2013b):
acos(1e7) = 0.0000 +16.8112i
acos(1e8) = 0.0000 +19.1138i

Octave 3.8.2:

acos(1e7) = 0.00000 - 16.81124i
acos(1e8) = 0.00000 - 19.11383i

Octave 4.0.0:

acos(1e7) = 0.00000 + 16.80543i
acos(1e8) = NaN + Infi

Conclusion: The changes from Octave 3.8.2 to 4.0.0

  • improved the Matlab compatibility regarding the sign of the imaginary part.
  • but partly destroyed the compatibility regarding the values of the real part. (for input values >~6e7)


(In my calculations I only used the real part, so that's why I bothered about this most...)

Anonymous
Thu 09 Jul 2015 11:08:25 PM UTC, comment #1: 

This is probably related to the correction to the phase mentioned in the NEWS file for 4.0.0


 ** The trigonometric functions asin and acos return different phase values
    from previous versions of Octave when the input is outside the principal
    branch ([-1, 1]).  If the real portion of the input is greater than 1 then
    the limit from below is taken.  If the real portion is less than 1 then the
    limit from above is taken.  This criteria is consistent with several other
    numerical analysis software packages.


This was done because 3.8.2 did not calculate the same result as Matlab.  Could you check your calculations in 3.8.2 because I don't believe they yield the same result as Matlab.  For 3.8.2 on Linux with the two calculations I get


acos(1e7)
ans =   0.00000 - 16.81124i
acos(1e8)
ans =   0.00000 - 19.11383i


where the sign of the imaginary portion is negative and the opposite of Matlab.


Rik <rik5>
Group administrator
Thu 09 Jul 2015 06:53:00 PM UTC, original submission:  

Here is a tiny "script" to show the behavior:


acos(1e7)
acos(1e8)


This gives you different results, depending where you run it (Octave 3.8.2, Octave 4.0.0, or Matlab):

  • The first calculation returns acos(1e7) = 0 + 16.8 i in all three programs.


  • The second calculation returnes acos(1e8) = 0 + 19.1 i in Matlab. It returnes the same value in Octave 3.8.2.


  • But in Octave 4.0.0 the result of the second line now is acos(1e8) = NaN + Inf i .


So Octave WAS behaving Matlab compatible here, but version 4.0.0 is NOT Matlab compatible any more in this case!

I don't know if acos is well defined for input values >1 at all. But Matlab decided to always returns something with a real part of 0, for big input numbers. Octave did the same, but now (version 4.0.0) it suprisingly "jumps" from a real part of 0 up to a real part of NaN somewhere between 1e7 and 1e8.

I have observed this change in behavior on the Linux release versions as well as on the Windows release versions (3.8.2-3 mxe inofficial built, 4.0.0 official release).

Hartmut <hardy>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #37692:  c99-acos-test.c added by mtmiller (307B - text/x-csrc)
file #37693:  cxx11-acos-test.cc added by mtmiller (351B - text/x-c++src)
file #37480:  complex-fns.diff added by mtmiller (6KiB - text/x-diff)
file #37476:  testx4.m added by caliari (2KiB - application/vnd.wolfram.mathematica.package)
file #37472:  complex-fns.diff added by mtmiller (6KiB - text/x-diff)
file #37293:  bug_45507.v4.patch added by lachlan (3KiB - text/x-diff)
file #36541:  bug45507.v3.diff added by mtmiller (3KiB - text/x-diff)
file #36198:  bug_45507.v2.cset added by lachlan (6KiB - application/octet-stream)
file #36167:  bug_45507.cset added by lachlan (2KiB - application/octet-stream)
file #34451:  testx3.m added by caliari (2KiB - application/vnd.wolfram.mathematica.package)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by lachlan (Updated the item)
  • -email is unavailable- added by caliari (Updated the item)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by hardy (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 24 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-09-02 rik5 StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2016-07-03 mtmiller Attached File- Added c99-acos-test.c, #37692
        Attached File- Added cxx11-acos-test.cc, #37693
    2016-07-03 lachlan StatusFixed In Progress
        Open/ClosedClosed Open
    2016-06-30 mtmiller StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2016-06-23 lachlan StatusFixed In Progress
    2016-06-23 lachlan Open/ClosedClosed Open
        Release4.0.0 dev
        Summaryacos returns different results on big input values acos returns different results on big input values on gcc 4.8
    2016-06-17 mtmiller StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2016-06-15 mtmiller Attached File- Added complex-fns.diff, #37480
    2016-06-14 caliari Attached File- Added testx4.m, #37476
    2016-06-14 mtmiller Attached File- Added complex-fns.diff, #37472
        StatusNeed Info In Progress
    2016-05-28 lachlan Attached File- Added bug_45507.v4.patch, #37293
    2016-03-04 mtmiller Attached File- Added bug45507.v3.diff, #36541
    2016-01-29 lachlan Attached File- Added bug_45507.v2.cset, #36198
    2016-01-26 lachlan Attached File- Added bug_45507.cset, #36167
    2015-07-16 caliari Attached File- Added testx3.m, #34451
    2015-07-09 rik5 StatusNone Need Info

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code