bugGNU Octave - Bugs: bug #49091, MinGW std::acosh less accurate...

 
 

bug #49091: MinGW std::acosh less accurate than Linux versions

Submitter:  Avinoam Kalma <avinoam>
Submitted:  Wed 14 Sep 2016 07:50:58 PM UTC
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Ready For Test Assigned to:  None
Originator Name:  Avinoam Open/Closed:  * Open
Release:  * dev Operating System:  * Microsoft Windows
* Mandatory Fields

Add a New Comment (Rich Markup)
   

Comment Type & Canned Response:
       

 

( Jump to the original submission )

Sun 07 Aug 2022 08:36:45 AM UTC, comment #28: 
Markus Mützel <mmuetzel>
Project Member
Fri 05 Aug 2022 05:37:14 PM UTC, comment #27: 

Turns out there were conditions that checked something like `x < 0.0`. But `-0.0 < 0.0` compares false, and the logic for the other side of the branch cut was used erroneously.

I pushed a patch to MXE Octave that will hopefully fix this:
https://hg.octave.org/mxe-octave/rev/9d1395334e68

I think this is the last issue in this report. Marking as ready for test.

Markus Mützel <mmuetzel>
Project Member
Fri 05 Aug 2022 06:27:27 AM UTC, comment #26: 

The negative zero seems to be the cause of the issue indeed. Thank you all for tracking this down. I didn't know that `-10i` in Octave actually meant `-(0,10)` as opposed to `(0,-10)`.
With a modified version of the test program in comment #15, I see:

$ ./test-acosh.exe
x = 0.0000000000000000e+00
y = 1.0000000000000000e+01
    cacosh(+x+1i*y) = 2.9982229502979698 + 1.5707963267948966i
    cacosh(+x-1i*y) = 2.9982229502979698 + -1.5707963267948966i
    cacosh(-x-1i*y) = 2.9982229502979760 + -1.5707963267948966i
std::acosh(+x+1i*y) = 2.9982229502979698 + 1.5707963267948966i
std::acosh(+x-1i*y) = 2.9982229502979698 + -1.5707963267948966i
std::acosh(-x-1i*y) = 2.998222950297976 + -1.5707963267948966i

That could probably be solved by using symmetries of the function also for the default implementation in MinGW's CRT (not only for large input values).

Markus Mützel <mmuetzel>
Project Member
Thu 04 Aug 2022 09:39:20 PM UTC, comment #25: 

Comment #24 does seem to indicate that the sign bit of the real part of the complex number has an effect even when the value of the real portion is zero.  That sure looks like a problem with the library.

One could test directly in the C++ code by doing

std::complex<double> zcpp (-x, -y)

which should create a complex value (-0, -10i).

Rik <rik5>
Project Administrator
Thu 04 Aug 2022 07:39:25 PM UTC, comment #24: 

>> format bit
>> 10i
ans = 0000000000000000000000000000000000000000000000000000000000000000  0100000000100100000000000000000000000000000000000000000000000000i
>> -10i
ans = 1000000000000000000000000000000000000000000000000000000000000000  1100000000100100000000000000000000000000000000000000000000000000i
>> complex(0,-10)
ans = 0000000000000000000000000000000000000000000000000000000000000000  1100000000100100000000000000000000000000000000000000000000000000i
>> -complex(0,10)
ans = 1000000000000000000000000000000000000000000000000000000000000000  1100000000100100000000000000000000000000000000000000000000000000i

>> real(acosh(10i))
ans = 0100000000000111111111000101110001010000011011010010101111011011
>> real(acosh(-10i))
ans = 0100000000000111111111000101110001010000011011010010101111101001
>> real(acosh(complex(0, -10)))
ans = 0100000000000111111111000101110001010000011011010010101111011011
>> real(acosh(-complex(0, 10)))
ans = 0100000000000111111111000101110001010000011011010010101111101001

Nicholas Jankowski <nrjank>
Project Member
Thu 04 Aug 2022 07:18:53 PM UTC, comment #23: 

Compare

format bit
-10i
complex (0, -10)
-complex (0, 10)
## etc.

Could having the sign bit set on the real part cause the problem here?

John W. Eaton <jwe>
Project Administrator
Thu 04 Aug 2022 06:35:45 PM UTC, comment #22: 

format bit; acosh(10i), acosh(-10i)
ans = 0100000000000111111111000101110001010000011011010010101111011011 0011111111111001001000011111101101010100010001000010110100011000i
ans = 0100000000000111111111000101110001010000011011010010101111101001 1011111111111001001000011111101101010100010001000010110100011000i

real parts do deviate

Nicholas Jankowski <nrjank>
Project Member
Thu 04 Aug 2022 06:29:25 PM UTC, comment #21: 

Could you try

$ octave --eval "format bit; acosh(10i), acosh(-10i)"

to verify that the computed values are different, not just the printed conversion from binary to decimal?

John W. Eaton <jwe>
Project Administrator
Thu 04 Aug 2022 05:17:01 PM UTC, comment #20: 

I tested with the following compiled to an .oct file:

#include <iostream>
#include <cmath>
#include <complex>

#include <octave/oct.h>
#include <octave/lo-specfun.h>

DEFUN_DLD(test_acosh_oct, args, nargout,
          "-*- texinfo -*-\n\
@deftypefn {} {@var{retval} =} test_acosh_oct (@var{input1}, @var{input2})\n\
@seealso{}\n\
@end deftypefn")
{
//  octave_value_list retval;
  int nargin = args.length ();
  double x = 0;
  double y = 10;
  std::cout.precision(17);
  std::cout << "x = " << x << std::endl;
  std::cout << "y = " << y << std::endl;

  // libstdc++
  std::complex<double> zcpp (x, y);
  std::complex<double> retcpp = std::acosh (zcpp);
  std::cout << "std::acosh(x+1i*y) = " << retcpp.real () << " + " << retcpp.imag () << "i" << std::endl;

  zcpp.imag (-y);
  retcpp = std::acosh (zcpp);
  std::cout << "std::acosh(x-1i*y) = " << retcpp.real () << " + " << retcpp.imag () << "i" << std::endl;

  // liboctave
  Complex zoct_p (x, y);
  Complex retoct_p = octave::math::acosh (zoct_p);
  Complex zoct_m (x, -y);
  Complex retoct_m = octave::math::acosh (zoct_m);

  // liboctinterp
  octave_value zoctint_p (zoct_p);
  octave_value zoctint_m (zoct_m);

  return ovl (retoct_p, retoct_m, zoctint_p.acosh (), zoctint_m.acosh ());
}

Executing it in Octave 7.2.0 yields:

>> mkoctfile test_acosh_oct.cc
>> [p, m] = test_acosh_oct;
x = 0
y = 10
std::acosh(x+1i*y) = 2.9982229502979698 + 1.5707963267948966i
std::acosh(x-1i*y) = 2.9982229502979698 + -1.5707963267948966i
>> [p, m, pi, mi] = test_acosh_oct;
x = 0
y = 10
std::acosh(x+1i*y) = 2.9982229502979698 + 1.5707963267948966i
std::acosh(x-1i*y) = 2.9982229502979698 + -1.5707963267948966i
>> fprintf ('acosh(x+iy) = %.16f + %.16fi\n', real (p), imag (p));
acosh(x+iy) = 2.9982229502979698 + 1.5707963267948966i
>> fprintf ('acosh(x-iy) = %.16f + %.16fi\n', real (m), imag (m));
acosh(x-iy) = 2.9982229502979698 + -1.5707963267948966i
>> fprintf ('acosh(x+iy) = %.16f + %.16fi\n', real (pi), imag (pi));
acosh(x+iy) = 2.9982229502979698 + 1.5707963267948966i
>> fprintf ('acosh(x-iy) = %.16f + %.16fi\n', real (mi), imag (mi));
acosh(x-iy) = 2.9982229502979698 + -1.5707963267948966i
>> z = acosh (10i);
>> printf ("Octave acosh(10i)  = %.16f + %.16fi\n", real (z), imag (z));
Octave acosh(10i)  = 2.9982229502979698 + 1.5707963267948966i
>> z = acosh (-10i);
>> printf ("Octave acosh(-10i) = %.16f + %.16fi\n", real (z), imag (z));
Octave acosh(-10i) = 2.9982229502979760 + -1.5707963267948966i

All results but for the last command look good to me.
I don't know what is going on here...

Is the parser doing something strange when it comes across `-10i`? 😖
The following result is fine:

>> z = acosh (complex(0,-10));
>> printf ("Octave acosh(-10i) = %.16f + %.16fi\n", real (z), imag (z));
Octave acosh(-10i) = 2.9982229502979698 + -1.5707963267948966i

Why only on Windows???

Markus Mützel <mmuetzel>
Project Member
Wed 03 Aug 2022 05:43:11 PM UTC, comment #19: 

Tracing the code, acosh is a mapper function in libinterp which forwards to liboctave.  In liboctave/numeric/lo-specfun.h there is a mapper of octave::math::acosh to

    inline Complex acosh (const Complex& x) { return std::acosh (x); }

This looks alright to me.

@Markus: What happens if you make a .oct file which has the C++ code that calls std::acosh and call that from Octave?

If that works, next step-wise change would be to use Complex and call octave::numeric::acosh directly.

Rik <rik5>
Project Administrator
Wed 03 Aug 2022 04:20:38 PM UTC, comment #18: 

That doesn't seem to make a difference:

octave:1> z = acosh (10i);
octave:2> printf ("Octave acosh(10i)  = %.16f + %.16fi\n", real (z), imag (z));
Octave acosh(10i)  = 2.9982229502979698 + 1.5707963267948966i
octave:3> z = acosh (-10i);
octave:4> printf ("Octave acosh(-10i) = %.16f + %.16fi\n", real (z), imag (z));
Octave acosh(-10i) = 2.9982229502979760 + -1.5707963267948966i

Markus Mützel <mmuetzel>
Project Member
Wed 03 Aug 2022 03:37:44 PM UTC, comment #17: 

it appears to be windows specific.  currently on both v6.4 and v7.2 I get:

>> z = acosh (10i);
>> printf ("Octave acosh(10i)  = %.16f + %.16fi\n", real (z), imag (z));
Octave acosh(10i)  = 2.9982229502979698 + 1.5707963267948966i
>> z = acosh (-10i);
>> printf ("Octave acosh(-10i) = %.16f + %.16fi\n", real (z), imag (z));
Octave acosh(-10i) = 2.9982229502979760 + -1.5707963267948966i

real difference is -6.2172e-15, or ~28eps.  while on v8 in linux the real difference is 'zero', matching your comment #16 result.

Nicholas Jankowski <nrjank>
Project Member
Wed 03 Aug 2022 03:01:28 PM UTC, comment #16: 

I think this is just rounding.  Instead of using "format long", use "printf" just as you would in C.

z = acosh (10i);
printf ("Octave acosh(10i)  = %.16f + %.16fi\n", real (z), imag (z));
z = acosh (-10i);
printf ("Octave acosh(-10i) = %.16f + %.16fi\n", real (z), imag (z));

On a Linux machine, the result for me is

Octave acosh(10i)  = 2.9982229502979698 + 1.5707963267948966i
Octave acosh(-10i) = 2.9982229502979698 + -1.5707963267948966i

Rik <rik5>
Project Administrator
Wed 03 Aug 2022 07:51:16 AM UTC, comment #15: 

I'm not sure what is happening for acosh(-10i) in Octave.

Test program:

#include <math.h>
#include <stdio.h>
#include <complex.h>

#include <iostream>
#include <cmath>
#include <complex>

int main()
{
  double x = 0;
  double y = 10;

  printf("x = %.16e\n", x);
  printf("y = %.16e\n", y);

  __complex__ double z;
  __real__ z = x;
  __imag__ z = y;

  // CRT
  __complex__ double ret = cacosh (z);
  printf("    cacosh(x+1i*y) = %.16f + %.16fi\n", __real__ ret, __imag__ ret);

  __imag__ z = -y;
  ret = cacosh (z);
  printf("    cacosh(x-1i*y) = %.16f + %.16fi\n", __real__ ret, __imag__ ret);

  // libstdc++
  std::complex<double> zcpp (x, y);
  std::complex<double> retcpp = std::acosh (zcpp);
  std::cout.precision(17);
  std::cout << "std::acosh(x+1i*y) = " << retcpp.real () << " + " << retcpp.imag () << "i" << std::endl;

  zcpp.imag (-y);
  retcpp = std::acosh (zcpp);
  std::cout << "std::acosh(x-1i*y) = " << retcpp.real () << " + " << retcpp.imag () << "i" << std::endl;

  return 0;
}

For MinGW, I see the following output:

Markus@Markus-PC MINGW64 ~
$ g++ test-acosh.cc -o test-acosh.exe

Markus@Markus-PC MINGW64 ~
$ ./test-acosh.exe
x = 0.0000000000000000e+00
y = 1.0000000000000000e+01
    cacosh(x+1i*y) = 2.9982229502979698 + 1.5707963267948966i
    cacosh(x-1i*y) = 2.9982229502979698 + -1.5707963267948966i
std::acosh(x+1i*y) = 2.9982229502979698 + 1.5707963267948966i
std::acosh(x-1i*y) = 2.9982229502979698 + -1.5707963267948966i

However, Octave still calculates the following:

$ octave --eval "format long; acosh(10i), acosh(-10i)"
ans =  2.998222950297970 + 1.570796326794897i
ans =  2.998222950297976 - 1.570796326794897i

Why does the simple C++ program calculate "good" values (using CRT or C++ functions)? But Octave still calculates "odd" values?

Markus Mützel <mmuetzel>
Project Member
Wed 03 Aug 2022 06:34:11 AM UTC, comment #14: 
Markus Mützel <mmuetzel>
Project Member
Wed 27 Jul 2022 08:02:55 AM UTC, comment #13: 
Markus Mützel <mmuetzel>
Project Member
Sat 23 Jul 2022 05:12:00 PM UTC, comment #12: 

I pushed a patch for the CRT of mingw-w64 to the default branch of MXE Octave here:
https://hg.octave.org/mxe-octave/rev/0825abaf61a7

With it, the tests with large inputs for inverse trigonometric functions pass for me with the same tolerances than on Linux (with libm from glibc?).

If that works for us, we could maybe send it upstream for review...

There are still failing tests for input that is pure negative imaginary.

Fwiw: Openlibm doesn't get that right either...

Markus Mützel <mmuetzel>
Project Member
Mon 09 Aug 2021 09:56:34 AM UTC, comment #11: 

These tests are still failing on Windows.
Afaict, the other bug was closed when the tests were changed to a xtests or marked as known issues with bug report numbers. IIUC, this report is tracking an actual fix. So, I think it should remain open.

Setting release to "dev".

Markus Mützel <mmuetzel>
Project Member
Mon 09 Aug 2021 09:13:28 AM UTC, comment #10: 

Per Rik this bug is suppose to be closed when bug #48312 is closed, this happened on 16-11-2016.

Anonymous
Thu 15 Sep 2016 09:05:28 PM UTC, comment #9: 

I don't think it matters too much.  As long as they are linked so that they get closed at the same time.

Rik <rik5>
Project Administrator
Thu 15 Sep 2016 06:45:55 PM UTC, comment #8: 

I suppose the acos FAILs mentioned in bug #48312 are the same.
Should this bug be listed there as dependency or the other way round?

Philip Nienhuis <philipnienhuis>
Project Member
Thu 15 Sep 2016 06:16:49 PM UTC, comment #7: 

I marked the failing acosh BIST tests with this bug number (http://hg.savannah.gnu.org/hgweb/octave/rev/4d9c371d2cca).  I also added a FIXME note saying that we are waiting for upstream to improve their library.  I've changed the bug title and marked it as postponed.

I'm also re-building with MXE to test that the failures are now reported correctly as known failures, but that will take a few hours.

Rik <rik5>
Project Administrator
Thu 15 Sep 2016 05:03:45 PM UTC, comment #6: 

I like your third option.  I don't feel we should have to dumb down Octave's tests for a single failing library on a single platform.  And I also don't like disabling a std:: library when there is a good possibility that it will be fixed upstream at some point.

Rik <rik5>
Project Administrator
Thu 15 Sep 2016 04:35:11 PM UTC, comment #5: 

A third option is to add this bug number to the tests, which reports them as XFAIL and this bug can remain open or be closed. Someone might want to report this to the MinGW project, if it's something that can be handled there.

I already marked some tests with bug #45507, which is about a similar failure on systems with an older gcc and/or libm system library. We are continuing to use the std::acosh function on these older GNU systems even though the results fail some tolerance tests.

Mike Miller <mtmiller>
Project Administrator
Thu 15 Sep 2016 04:06:22 PM UTC, comment #4: 

The results I am seeing in fntests.log are

***** test
 re = 2.99822295029797;
 im = pi/2;
 assert (acosh (-10i), re - i*im);
!!!!! test failed
ASSERT errors for:  assert (acosh (-10i),re - i * im)

  Location  |  Observed  |  Expected  |  Reason
     ()      2.9982-1.5708i 2.9982-1.5708i   Abs err 6.2172e-015 exceeds tol 0
shared variables
  scalar structure containing the fields:

    rt2 =  1.4142
    rt3 =  1.7321
***** test
 re = single (2.99822295029797);
 im = single (pi/2);
 assert (acosh (single (10i)), re + i*im, 5*eps ("single"));
 assert (acosh (single (-10i)), re - i*im, 5*eps ("single"));
!!!!! test failed
ASSERT errors for:  assert (acosh (single (-10i)),re - i * im,5 * eps ("single"))

  Location  |  Observed  |  Expected  |  Reason
     ()      2.9982-1.5708i 2.9982-1.5708i   Abs err 7.1526e-006 exceeds tol 5.9605e-007
shared variables
  scalar structure containing the fields:

    rt2 =  1.4142
    rt3 =  1.7321

When using doubles, it is failing by 31 eps.  When using floats it is failing by60 eps ("single").  Those are pretty big tolerances.  Do we want to add those in or disable the library on Windows?

Rik <rik5>
Project Administrator
Thu 15 Sep 2016 12:26:11 AM UTC, comment #3: 

I checked and _have_feature_ reports that I am using the std library version.  So, that kind of sucks.  The std library is less accurate than direct computation on Windows machines.

Do we want to change the library or the test?  We could certainly add to the #ifdef.

#if defined (HAVE_COMPLEX_STD_ACOSH) && ! defined (OCTAVE_USE_WINDOWS_API)
      return std::acosh (x);
#else
      return log (x + sqrt (x + 1.0) * sqrt (x - 1.0));
#endif

The second option is to increase the tolerance of the test.

Rik <rik5>
Project Administrator
Wed 14 Sep 2016 10:23:48 PM UTC, comment #2: 

The mxe-octave binary distribution should include config.log in the share/octave/4.2.0-rc1/etc directory, that should tell you whether std::acosh was found by configure. Also running

>> __have_feature__ COMPLEX_STD_ACOSH
ans = 1

in the interpreter should tell the same thing.

Mike Miller <mtmiller>
Project Administrator
Wed 14 Sep 2016 10:17:54 PM UTC, comment #1: 

I guess the question is why acosh (10i) is different from acosh (-10i).  On Linux they are the same.  But on Windows, I get

>> acosh (10i)
ans =  2.99822295029797 + 1.57079632679490i
>> acosh (-10i)
ans =  2.99822295029798 - 1.57079632679490i

In liboctave/libnumeric/lo-specfun.cc the acosh function is defined as

    Complex
    acosh (const Complex& x)
    {
#if defined (HAVE_COMPLEX_STD_ACOSH)
      return std::acosh (x);
#else
      return log (x + sqrt (x + 1.0) * sqrt (x - 1.0));
#endif
    }

If I execute the #else path in Octave I get the right answer.

x =   0 + 10i
>> log (x + sqrt (x + 1.0) * sqrt (x - 1.0))
ans =  2.99822295029797 + 1.57079632679490i
>> x = -10i
x =  -0 - 10i
>> log (x + sqrt (x + 1.0) * sqrt (x - 1.0))
ans =  2.99822295029797 - 1.57079632679490i

I don't know how MXE Octave was build, but this suggests that the std library might be at fault on Windows.

Rik <rik5>
Project Administrator
Wed 14 Sep 2016 07:50:58 PM UTC, original submission:  

>> test mappers.cc-tst
***** xtest
 x = [1, -1, i, -i] .* 1e150;
 v = [0, pi, pi/2, pi/2];
 assert (real (acos (x)), v);
!!!!! known failure
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
***** test
 re = 2.99822295029797;
 im = pi/2;
 assert (acosh (-10i), re - i*im);
!!!!! test failed
ASSERT errors for:  assert (acosh (-10i),re - i * im)

  Location  |  Observed  |  Expected  |  Reason
     ()      2.9982-1.5708i 2.9982-1.5708i   Abs err 6.2172e-015 exceeds tol 0
shared variables
  scalar structure containing the fields:

    rt2 =  1.4142
    rt3 =  1.7321

Using a tolerance of 7e-15 instead of zero will solve this.

Avinoam Kalma <avinoam>
Project Member

 

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

Attach Files:
   
   
Comment:
   

No files currently attached

 

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 nrjank (Posted a comment)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by avinoam (Submitted the item)
  • -email is unavailable- added by avinoam
  •  

    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 project members can vote.

     

     

     

     

    Follow 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-08-05 mmuetzel StatusPostponed Ready For Test
    2021-08-09 mmuetzel Release4.2.0 dev
    2016-11-17 mtmiller Release4.2.0-rc1 4.2.0
    2016-09-15 rik5 StatusConfirmed Postponed
        SummaryMXE Octave 4.20-rc1: test mappers.cc-tst fails MinGW std::acosh less accurate than Linux versions
    2016-09-14 rik5 StatusNone Confirmed
    2016-09-14 avinoam Carbon-Copy- Added -email is unavailable-

    Back to the top


    Powered by Savane 3.9