bugGNU Octave - Bugs: bug #44310, log and functions which use log...

 
 

bug #44310: log and functions which use log return incorrect result

Submitter:  Michael Godfrey <godfrey>
Submitted:  Thu 19 Feb 2015 02:28:11 PM UTC
   
 
Category:  Libraries Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  godfrey Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Fri 17 Jun 2016 09:54:46 PM UTC, comment #27: 

This bug is now fixed on the default branch. Thanks for your patience and feedback. I've added some minimal unit testing of atanh to attempt to confirm correct output when given a large magnitude argument.

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

Mike Miller <mtmiller>
Group Member
Sun 06 Mar 2016 02:28:22 AM UTC, comment #26: 

Mike:  I tried it and it works as expected. It would be
great if the same fix could be applied to all the affected
functions (i.e. all that use the (unstable) functions of
log).

This should also improve preformance, but I have not
checked this. Anyhow, the key result is getting the correct
answer.

Thanks VERY much for doing this essentially as soon as
the C++11 was enabled!!

Michael Godfrey <godfrey>
Group Member
Fri 04 Mar 2016 12:00:40 AM UTC, comment #25: 

Please try the attached patch on the default branch. Now that we have C++11 enabled we can test for and use these functions if available.

If this patch works for you I will expand it to also address bug #45507 similarly.

(file #36540)

Mike Miller <mtmiller>
Group Member
Wed 27 Jan 2016 02:35:52 AM UTC, comment #24: 

This bug has been outstanding for quite a while.
The patch provided fixes the specific problem
with atanh, but it would be VERY good if someone
with the required C++ skills could implement the
suggested fix for the other elementary functions
that currently use numerically unstable transformations.

Michael Godfrey <godfrey>
Group Member
Thu 02 Jul 2015 01:27:33 AM UTC, comment #23: 

Attached is apatch which just fixes the reported
problem with atanh(0 -10i) and provides a test
for this case.
the elementary functions really should be reviewed
and upgraded, but this seems to depend on C++ and
possibly other issues.
this report can be closed, but another one about
the review should e opened.



(file #34358)

Michael Godfrey <godfrey>
Group Member
Thu 28 May 2015 04:21:51 PM UTC, comment #22: 

Attached is a changeset which adds an xtest
to mappers.cc.

This xtest checks for correct result of catanh.

Applying this patch willsimplyeo remind people
that this error requires fixing someday.



(file #34117)

Michael Godfrey <godfrey>
Group Member
Thu 30 Apr 2015 12:41:25 PM UTC, comment #21: 


As I mentioned in bug #41061 the current devel system
compiles correctly (on Fedora 21) using C++11, but with
a number of warnings. The only patch needed is simply change line
86 of stddef.h from #if !0 to #if 0. This will have to be made
conditional if there really are still some systems without a
definition of max_align_t.

Also, there are a number of other elementary functions in
mappers.cc which use the formal definition in terms of log,
but which have specific numerical codes in libm. These should
also be changed.

Michael Godfrey <godfrey>
Group Member
Thu 30 Apr 2015 09:26:20 AM UTC, comment #20: 

Mike: are you sure that atanh (in libm) requires C++11? Other
functions, such as log, are currently used in mappers.cc.
In any case, using the libm atanh is the right idea. And, as
llvm 3.5 needs C++11, conditionally using C++11 will, it appears,
be necessary at some point. is there a better (simpler) choice?

Michael Godfrey <godfrey>
Group Member
Thu 30 Apr 2015 02:31:08 AM UTC, comment #19: 

Adding dependency on bug #44965. I think the cleanest solution to this bug is to conditionally use the std::atanh() C++ function template, which is only available in C++11. If the compiler does not support C++11 or does not have std::atanh, fall back to the current method.

Mike Miller <mtmiller>
Group Member
Fri 20 Mar 2015 03:36:59 PM UTC, comment #18: 

Rik:  Would you consider adding the test which is in the
attached file #33194 (it should be made an XFAIL I suppose).

This will remind people that this error in atanh needs
to be fixed someday.



Michael Godfrey <godfrey>
Group Member
Tue 10 Mar 2015 02:47:09 PM UTC, comment #17: 

OK. I agree. Wrong answers do make me nervous!
But, 4.0 release should take precedence.

And, as you mention, the effect on users will vary depending
on what glibc they have.  And therefore, of course, it
affects people using other software that relies on
the log definition of actanh and other functions.

So, OK, but try to remember for sometime after 4.0 release.
At that point it might make sense to review the elementary
functions in mappers.cc as well. And, does Octave support
any other math library on systems which do not have glibc?
I did not think it did, but have not checked.

Michael Godfrey <godfrey>
Group Member
Tue 10 Mar 2015 02:17:34 PM UTC, comment #16: 

I don't agree with the characterization of this as a release blocker. It's a degradation in numerical accuracy to be sure, that may or may not be present in Octave depending on your system libraries. If I compare Octave 3.4 and Octave 4.0, I'll see the same result, so it's not a degradation introduced by the new release.

In addition, the fix may work for some users but may break the build for others. Yes, on the surface the change is as simple as you suggest.

However, the catanh and corresponding std::atanh(std::complex) functions are not generally available on all systems, so I think some thought should be put into how they are detected and conditionally used and what compiler flags we might have to introduce to make it work for all users.

Mike Miller <mtmiller>
Group Member
Sun 08 Mar 2015 09:41:03 PM UTC, comment #15: 

Mike:  I thought about this a bit more: How much work
is it for a C++ person to fix this by use of replacing

return log ((1.0 + x) / (1.0 - x)) / 2.0;
with
return catanh( x );

If it is not too much work to get all the declarations
straight, it would be good to get this into 4.0.

It is the following:
1. Regression
2. Incorrect result
3. Incompatible with Matlab

these formally qualify it as a blocker.
I f I could help more I would. If you tell me
what to do, I will test it.

Michael Godfrey <godfrey>
Group Member
Sun 08 Mar 2015 03:49:44 AM UTC, comment #14: 

Mike: You obviously know a lot about this.
I know lots of programming languages but C++ is not,
and never will be, one of them.

So, could you keep in mind that sometime after the 4.0 release
you will deal with this?

It does seem to me that catanh should use the glibc library,
as should other of the standard elementary functions.

Relying on the use of log is not, numerically, a good idea.

Michael Godfrey <godfrey>
Group Member
Sun 08 Mar 2015 01:22:11 AM UTC, comment #13: 

I cannot spell C++, so this is as far as I got:

numeric/lo-mappers.cc: In function 'Complex atanh(const Complex&)':
numeric/lo-mappers.cc:242:11: error: 'catanh' is not a member of 'std'
   return  std::catanh( x );
           ^
numeric/lo-mappers.cc:242:11: note: suggested alternative:
In file included from numeric/lo-mappers.cc:32:0:
numeric/lo-mappers.h:175:32: note:   'catanh'
 extern OCTAVE_API FloatComplex catanh (const FloatComplex& x);
                                ^
numeric/lo-mappers.cc:243:1: warning: control reaches end of non-void function [-Wreturn-type]
 }
 ^

Michael Godfrey <godfrey>
Group Member
Sun 08 Mar 2015 12:02:41 AM UTC, comment #12: 

Since Octave uses C++ wherever possible, you should look at the C++ function template std::atanh instead of catanh/catanhf. This is via the standard complex header file, so no need to include anything more.

This is also a relatively new feature of C++11, so there would need to be some Autoconf work done to make sure the function is available and fall back to the current equivalent log expression when the compiler doesn't provide the function.

Mike Miller <mtmiller>
Group Member
Sat 07 Mar 2015 10:56:30 PM UTC, comment #11: 

Where is the declaration for catanh?
Or, does it need an Octave declaration of some sort?
And, a hint about the declaration would help.

log seems to work without any additional declaration...

Michael Godfrey <godfrey>
Group Member
Sat 07 Mar 2015 09:19:10 PM UTC, comment #10: 

Good.

Would replacing 
  return log ((1.0 + x) / (1.0 - x)) / 2.0;
with
  return catanh( x );

just work?


Michael Godfrey <godfrey>
Group Member
Sat 07 Mar 2015 08:50:18 PM UTC, comment #9: 
Mike Miller <mtmiller>
Group Member
Sat 07 Mar 2015 08:09:52 PM UTC, comment #8: 

I tried looking at mappers.cc and its header .h files,
but I could not follow its logic for processing
atanh() for complex arguments. However, it appears to
me that it would be better in any case if atanh called the
appropriate libm function (atanh() or catanh()).

Making use the the log definition is mathematically correct
but may not be the best choice for numerical accuracy.
Obviously, the glibc folks think that catanh( ) should
be evaluated differently from using clog(xp/xm)/2.0.

The fact that long double complex anl = clogl(xp/xm)/2.0;
has "always" returned 4.672911e-17-1.471128e+00i suggests
that they may consider also returning this as the answer
for double complex anl = clogl(xp/xm)/2.0 to be more accurate
and this is likely true.  So, they are not likely to agree
to a change
in their code.

So, could a mappers.cc expert take a look at this, please?


Michael Godfrey <godfrey>
Group Member
Thu 26 Feb 2015 07:47:31 PM UTC, comment #7: 

Note that this is also a Matlab incompatbility.
Matlab gets the right answer.

Michael Godfrey <godfrey>
Group Member
Thu 26 Feb 2015 02:52:43 PM UTC, comment #6: 

Here is the conclusion of where the atanh error comes
from, and a patch is attached that tests for the error
in mappers.cc.


On FC14 Octave 3.8.1
octave:10> xp/xm
ans = -0.980198019801980 + 0.198019801980198i
octave:11> log(xp/xm)/2
ans =  0.000000000000000 + 1.471127674303735i

On FC14 Octave 3.9.1+
octave:4> xp/xm
ans = -0.980198019801980 + 0.198019801980198i
octave:5> log(xp/xm)/2
ans =  0.000000000000000 + 1.471127674303735i

On FC21 Octave 3.9.1+
octave:14> xp/xm
ans = -0.980198019801980 + 0.198019801980198i
octave:15> log(xp/xm)/2
ans =  4.67173054916526e-17 + 1.47112767430373e+00i


So, the error is in the log function in libm-2.20.
The error was introduced after libm-2.18.

My attempt to communicate this error to the glibc folks
appears not to have gotten through.

Applying the attached patch will at least remind everyone
that the error exists.


(file #33194)

Michael Godfrey <godfrey>
Group Member
Wed 25 Feb 2015 10:40:46 PM UTC, comment #5: 

I left out what I meant to include "below"
That is:


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

int main(void)
{
    float complex z = catanh(0-10i);
    printf("catanh(0-10i) = %e%+ei\n", creal(z), cimag(z));

    double complex x  = (0-10i);
    double complex xm = (1.0 - x);
    double complex xp = (1.0 + x);
    double complex an = clog(xp/xm)/2.0;
    long double complex anl = clogl(xp/xm)/2.0;
    double complex ans   = (clog(xp) - clog(xm))/2.0;
    printf("catanh (an)  =  %e%+ei\n", creal(an), cimag(an));
    printf("catanh (anl) =  %e%+ei\n", creal(anl), cimag(anl));
    printf("catanh (ans) =  %e%+ei\n", creal(ans), cimag(ans));
}
which gives:
catanh(0-10i) = 0.000000e+00-1.471128e+00i
catanh (an)  =  4.671731e-17-1.471128e+00i
catanh (anl) =  4.671731e-17-1.471128e+00i
catanh (ans) =  0.000000e+00-1.471128e+00i

This shows that clog(xp) - clog(xm) yields a different
and correct answer instead of clog(xp/xm)

Michael Godfrey <godfrey>
Group Member
Wed 25 Feb 2015 10:34:26 PM UTC, comment #4: 

The problem appears to be with the line:
 double complex an = clog(xp/xm)/2.0;
replacing this by:
 double complex ans   = (clog(xp) - clog(xm))/2.0;
corrects the problem with atanh(0 -10i)
as shown below. Using 2 calls to clog is less
efficient, but gets the right answer. I do not
know whether the problem came about due to a change in
the code produced for for xp/xm, or replacing clog(xp) - clog(xm)
by clog(xp/xm).

The attached patch tests for this error in mappers.cc


(file #33188)

Michael Godfrey <godfrey>
Group Member
Thu 19 Feb 2015 07:01:54 PM UTC, comment #3: 

@Michael: yes, it was just to be sure.

On Mageia-4 (also glibc):

>> format long
>> atanh(0 -10i)
ans =  4.67173054916526e-17 - 1.47112767430373e+00i
>>


Philip Nienhuis <philipnienhuis>
Group Member
Thu 19 Feb 2015 03:40:49 PM UTC, comment #2: 

I would not have expected this to affect Windows
since it is a glibc-2.20 problem (glibc provides
libm-2.20.so which is where log function exists).
It affects Fedora systems and probably other linux
based systems.

Michael Godfrey <godfrey>
Group Member
Thu 19 Feb 2015 02:56:01 PM UTC, comment #1: 

FTR, on Windows, 3.9.1+ from Feb 18 (w. 64-bit indexing):

>> atanh(0 -10i)
ans =  0.00000 - 1.47113i
>> format long
>> atanh(0 -10i)
ans =  0.000000000000000 - 1.471127674303735i
>>


Philip Nienhuis <philipnienhuis>
Group Member
Thu 19 Feb 2015 02:28:11 PM UTC, original submission:  

Example:

octave:1> atanh(0 -10i)
ans =  4.6717e-17 - 1.4711e+00i
octave:2>

This occurs on Fedora systems, and maybe others.
The library: /usr/lib64/libm-2.20.so contains the log function. Until fairly
recently this seemed to work correctly. But in  libm-2.20.so the following
happens using:


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

int main(void)
{
    float complex z = catanh(0-10i);
    printf("catanh(0-10i) = %e%+ei\n", creal(z), cimag(z));

    double complex x  = (0-10i);
    double complex xm = (1.0 - x);
    double complex xp = (1.0 + x);
    double complex an = clog(xp/xm)/2.0;
    long double complex anl = clogl(xp/xm)/2.0;
    printf("catanh (an)  =  %e%+ei\n", creal(an), cimag(an));
    printf("catanh (anl) =  %e%+ei\n", creal(an), cimag(an));
}
gives:
[godfrey@pbdsl4 atanh]$ gcc -lm -o testt atanh_t.cc
[godfrey@pbdsl4 atanh]$ testt
catanh(0-10i) = 0.000000e+00-1.471128e+00i
catanh (an)  =  4.671731e-17-1.471128e+00i
catanh (anl) =  4.671731e-17-1.471128e+00i
==================

Note that both clog and clogl fail.

Under libm-2.13.so the result for clog was:
 
catanh(0-10i) = 0.000000e+00-1.471128e+00i
catanh =  0.000000e+00-1.471128e+00i

While  4.671731e-17 is a small number (less than a usual eps) it
is not correct. This result causes the current Octave system to produce
incorrect results as it uses clog in its math function library.

This first showed up when using the current Fedora 21 x86_64 system.

I have attempted to report this to
 sourceware.org/bugzilla but I am not sure that the report
got through. In any case the error causes errors in Octave.

Michael Godfrey <godfrey>
Group Member

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #36540:  bug44310.diff added by mtmiller (2KiB - text/x-diff)
file #34358:  fix_atanh.diff added by godfrey (1KiB - text/x-patch)
file #34117:  mapper.diff added by godfrey (786B - text/x-patch)
file #33194:  complex_atanh_error.diff added by godfrey (678B - text/x-patch)
file #33188:  complex_atanh_error.diff added by godfrey (678B - text/x-patch)

 

Digest:
   bug dependencies.

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by godfrey (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 12 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2016-06-17 mtmiller StatusPatch Submitted Fixed
        Open/ClosedOpen Closed
    2016-03-04 mtmiller Attached File- Added bug44310.diff, #36540
        StatusConfirmed Patch Submitted
    2015-07-02 godfrey Attached File- Added fix_atanh.diff, #34358
    2015-05-28 godfrey Attached File- Added mapper.diff, #34117
    2015-04-30 mtmiller Dependencies- Depends on bugs #44965
    2015-03-07 mtmiller CategoryOctave Function Libraries
        Item GroupRegression Inaccurate Result
    2015-02-26 godfrey Attached File- Added complex_atanh_error.diff, #33194
    2015-02-25 godfrey Attached File- Added complex_atanh_error.diff, #33188
    2015-02-19 philipnienhuis StatusNone Confirmed

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code