bugGNU Octave - Bugs: bug #42583, log2() returns inaccurate result...

 
 

bug #42583: log2() returns inaccurate result for many integer powers of 2, unlike Matlab

Submitter:  Falk Tannhäuser <fata>
Submitted:  Wed 18 Jun 2014 11:59:41 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  Falk Tannhäuser Open/Closed:  * Closed
Release:  * 3.8.1 Operating System:  * Microsoft Windows
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 27 Feb 2023 02:39:34 AM UTC, comment #27: 

I've had bad experiences with the '-ffast-math' option.  The documentation from gcc for -ffast-math is:


       -ffast-math
           Sets the options -fno-math-errno, -funsafe-math-optimizations, -ffinite-math-only,
           -fno-rounding-math, -fno-signaling-nans, -fcx-limited-range and -fexcess-precision=fast.

           This option causes the preprocessor macro "__FAST_MATH__" to be defined.

           This option is not turned on by any -O option besides -Ofast since it can result in incorrect
           output for programs that depend on an exact implementation of IEEE or ISO
           rules/specifications for math functions. It may, however, yield faster code for programs that
           do not require the guarantees of these specifications.


It seems possible that having even one module compiled with -ffast_math might contaminate others.

Rik <rik5>
Group administrator
Sun 26 Feb 2023 03:00:12 PM UTC, comment #26: 

I'm not sure this is related, but some years ago I was having some strange crashes and it turned out that the issue was related
 to linking to two different versions of BLAS and LAPACK at the same time because the qrupdate I was linking to was statically linked  to a different version than that I was explicitly linking to ...
Could there be a similar problem here?

Carlo de Falco <cdf>
Group Member
Sun 26 Feb 2023 12:31:03 AM UTC, comment #25: 

I thought it is lto, but I recompiled with default flags and it
was still a problem. Anyway, the qrupdate was compiled with "-O3 --ffast-math" (those are the original flags). I recompiled with -O2 and now it looks fine.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Sun 26 Feb 2023 12:23:45 AM UTC, comment #24: 

I don't understand either.  The function std::exp2 shouldn't have any relation to qrupdate library.

Have you tried enabling 64-bit libraries, but no other optimizations besides -O2?  My guess would still be that one of the tuning options such as '-march=native' or '-flto' might cause this.

Rik <rik5>
Group administrator
Sat 25 Feb 2023 10:56:55 PM UTC, comment #23: 

It looks like qrupdate causes the trouble. If I disable it,
64-bit works fine. I do not understand how it can be happening.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Sat 25 Feb 2023 10:45:36 PM UTC, comment #22: 

I recompiled w/o all those 64-bit Fortran libraries and now
64-bit version is correct.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Sat 25 Feb 2023 08:58:48 PM UTC, comment #21: 

I compiled with default optimization flags both 64-bit pointer and 32-bit pointer. The result is the same. With 32-bit pointers I got the same answer as Rik:

octave:1> __octave_config_info__ ("config_opts")
ans =
octave:2> pow2 (-1026:-1020)
ans =

  1.4e-309  2.8e-309  5.6e-309  1.1e-308  2.2e-308  4.5e-308  8.9e-308

With 64-bit pointers:

ctave:1> __octave_config_info__ ("config_opts")
ans =  '--enable-64' '--with-blas=flexiblas64' '--with-qrupdate=qrupdate64' '--with-suitesparseconfig=suitesparseconfig64' '--with-amd=amd64' '--with-camd=camd64' '--with-colamd=colamd64' '--with-ccolamd=ccolamd64' '--with-cholmod=cholmod64' '--with-spqr=spqr64' '--with-cxsparse=cxsparse64' '--with-umfpack=umfpack64' '--with-klu=klu64' '--with-arpack=arpack64' '--without-sundials_nvecserial' '--without-sundials_ida' '--without-sundials_sunlinsolklu'
octave:2> pow2 (-1026:-1020)
ans =

   0   0   0   0  2.2e-308  4.5e-308  8.9e-308


I also tried clang (32-bit and 64-bit) -- the same discrepancy.

Trying smaller steps:


octave:37> pow2 (-(1022+257*eps))
ans = 0
octave:38> pow2 (-(1022+256*eps))
ans = 2.2251e-308
octave:39> format bit
octave:40> -(1022+257*eps)
ans = 1100000010001111111100000000000000000000000000000000000000000001
octave:41> -(1022+256*eps)
ans = 1100000010001111111100000000000000000000000000000000000000000000


Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Sat 25 Feb 2023 08:37:21 PM UTC, comment #20: 

The code for pow2 is in libinterp/corefcn/pow2.cc.  For doubles,


for (octave_idx_type i = 0; i < y.numel (); i++)
  y.xelem (i) = std::exp2 (x.xelem (i));


So, Octave is just making use of the C++ math library.

On my machine,


pow2 (-1026:-1020)
ans =

  1.4e-309  2.8e-309  5.6e-309  1.1e-308  2.2e-308  4.5e-308  8.9e-308


But this is just with very simple compile flags "-O2 -pipe".


Rik <rik5>
Group administrator
Sat 25 Feb 2023 05:21:10 PM UTC, comment #19: 

The test

test libinterp/corefcn/data.cc-tst
***** assert <*42583> (all (log2 (pow2 (-1074:1023)) == -1074:1023))
!!!!! regression: https://octave.org/testfailure/?42583
assert (all (log2 (pow2 (-1074:1023)) == -1074:1023)) failed


is failing on octave 9 (c80cf1588ed0) compiled with 64-bit pointers:


octave:12> log2 (pow2 (-1026:-1020))
ans =

   -Inf   -Inf   -Inf   -Inf  -1022  -1021  -1020

octave:13> pow2 (-1026:-1020)
ans =

   0   0   0   0  2.2e-308  4.5e-308  8.9e-308

octave:14> format bit; pow2 (-1026:-1020)
ans =

 Columns 1 and 2:

  0000000000000000000000000000000000000000000000000000000000000000  0000000000000000000000000000000000000000000000000000000000000000

 Columns 3 and 4:

  0000000000000000000000000000000000000000000000000000000000000000  0000000000000000000000000000000000000000000000000000000000000000

 Columns 5 and 6:

  0000000000010000000000000000000000000000000000000000000000000000  0000000000100000000000000000000000000000000000000000000000000000

 Column 7:

  0000000000110000000000000000000000000000000000000000000000000000


Default compile (same computer, gcc, and compiler flags) passes.


> __octave_config_info__ ("config_opts")
ans =  '--enable-64' '--with-blas=flexiblas64' '--with-qrupdate=qrupdate64' '--with-suitesparseconfig=suitesparseconfig64' '--with-amd=amd64' '--with-camd=camd64' '--with-colamd=colamd64' '--with-ccolamd=ccolamd64' '--with-cholmod=cholmod64' '--with-spqr=spqr64' '--with-cxsparse=cxsparse64' '--with-umfpack=umfpack64' '--with-klu=klu64' '--with-arpack=arpack64' '--without-sundials_nvecserial' '--without-sundials_ida' '--without-sundials_sunlinsolklu' 'CFLAGS=-ggdb3 -O2 -march=native -mtune=native -mavx -mavx2 -flto=auto' 'CXXFLAGS=-ggdb3 -O2 -march=native -mtune=native -mavx -mavx2 -flto=auto' 'FFLAGS=-ggdb3 -O2 -march=native -mtune=native -mavx -mavx2 -flto=auto' 'PKG_CONFIG_PATH=/opt/rh/gcc-toolset-12/root/usr/lib64/pkgconfig'


Tried with both gcc version 11.3.1 20221121 (Red Hat 11.3.1-4) (GCC), and gcc version 12.2.1 20221121 (Red Hat 12.2.1-7) (GCC).

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Sun 22 Jun 2014 12:59:30 AM UTC, comment #18: 

Great.  I'm marking this bug as fixing and closing the report.

Rik <rik5>
Group administrator
Sat 21 Jun 2014 07:25:52 PM UTC, comment #17: 

@Rik (in response to https://savannah.gnu.org/bugs/?42583#comment10)
I succeeded in downloading all required tools and libraries (except Java) and in compiling Octave under Cygwin (64 bit) from the Mercurial sources, using GCC 4.9.0. The problem is gone away both in the stable branch (version 3.8.2-rc1) and in the default branch (version 4.1.0+).
"make check" in the stable branch gives the following summary:


  PASS     11547
  FAIL        21
  XFAIL        7
  SKIPPED     43


According to test/fntests.log, concerning log2() there only subsist assert failures for infinite complex arguments, a case probably rarely encountered in practice.

Falk Tannhäuser <fata>
Fri 20 Jun 2014 11:04:49 AM UTC, comment #16: 

@Stefan Mahr: yes, sorry. Bootstrapping solved my problem.

Kai Torben Ohlhus <siko1056>
Group Member
Fri 20 Jun 2014 08:47:49 AM UTC, comment #15: 

I just googled a bit and found this macro to provide gnulib functions in their namespace:

#define GNULIB_NAMESPACE gnulib

https://www.gnu.org/software/gnulib/manual/html_node/A-C_002b_002b-namespace-for-gnulib.html#A-C_002b_002b-namespace-for-gnulib

Kai Torben Ohlhus <siko1056>
Group Member
Fri 20 Jun 2014 08:46:55 AM UTC, comment #14: 

You most probably forgot to ./bootstrap

Stefan Mahr <dac922>
Fri 20 Jun 2014 08:40:23 AM UTC, comment #13: 

The past changes resulted in a compilation error ('log2' is not a member of 'gnulib') for my Debian Jessie with GCC 4.9:

../liboctave/numeric/lo-mappers.cc:95:10: error: 'log2' is not a member of 'gnulib'
   return gnulib::log2 (x);
          ^
../liboctave/numeric/lo-mappers.cc:95:10: note: suggested alternative:
In file included from /usr/include/features.h:374:0,
                 from /usr/include/stdint.h:25,
                 from /usr/lib/gcc/x86_64-linux-gnu/4.9/include/stdint.h:9,
                 from ../oct-conf-post.h:167,
                 from ../config.h:3372,
                 from /
../liboctave/numeric/lo-mappers.cc:25:
/usr/include/x86_64-linux-gnu/bits/mathcalls.h:144:1: note:   'log2'
 __MATHCALL (log2,, (Mdouble __x));
 ^
../liboctave/numeric/lo-mappers.cc: In function 'float xlog2(float)':
../liboctave/numeric/lo-mappers.cc:310:10: error: 'log2f' is not a member of 'gnulib'
   return gnulib::log2f (x);
          ^
../liboctave/numeric/lo-mappers.cc:310:10: note: suggested alternative:
In file included from ../libgnu/math.h:27:0,
                 from /usr/include/c++/4.9/cmath:44,
                 from /usr/include/c++/4.9/complex:44,
                 from ../liboctave/util/oct-cmplx.h:27,
                 from ../liboctave/numeric/lo-mappers.h:29,
                 from ../liboctave/numeric/lo-mappers.cc:32:
/usr/include/x86_64-linux-gnu/bits/mathcalls.h:144:1: note:   'log2f'
 __MATHCALL (log2,, (Mdouble __x));
 ^

Kai Torben Ohlhus <siko1056>
Group Member
Fri 20 Jun 2014 12:10:02 AM UTC, comment #12: 

I also added a test for this problem:

http://hg.savannah.gnu.org/hgweb/octave/rev/9d185537e5d1

John W. Eaton <jwe>
Group administrator
Thu 19 Jun 2014 11:35:51 PM UTC, comment #11: 

I sent a mail to the Cygwin list: https://cygwin.com/ml/cygwin/2014-06/msg00314.html
Although I found some discussions about log2 in the list archives, I didn't see anybody mentioning accuracy issues.
BTW I got Octave for Cygwin through the usual setup-x86_64.exe installation, so I suppose it was native compilation.

Falk Tannhäuser <fata>
Thu 19 Jun 2014 11:29:09 PM UTC, comment #10: 

@Falk: jwe added the gnulib module for log2 to Octave (http://hg.savannah.gnu.org/hgweb/octave/rev/ff4da3c8ed16).  Can you try compiling from Mercurial source code and see if the problem has gone away?  The change was made on the stable branch which will become 3.8.2 fairly soon.

Rik <rik5>
Group administrator
Thu 19 Jun 2014 05:53:32 PM UTC, comment #9: 

We don't cross compile for Cygwin, at least that I know of.  It's not something that the mxe-octave build system does.


John W. Eaton <jwe>
Group administrator
Thu 19 Jun 2014 05:52:36 PM UTC, comment #8: 

Hmm, I just looked at and ran the bootstrap script and it looks like we aren't using the gnulib log2 or log2-ieee module.  Since the latter depends on log2, I think it is best to use the log2-ieee module.  I'm looking at making that change now.

Reporting the problem for the Cygwin library would be OK, but since this is known due to the fixes in gnulib, I suspect it would just be a duplicate report.

John W. Eaton <jwe>
Group administrator
Thu 19 Jun 2014 05:38:37 PM UTC, comment #7: 

That is probably the best thing to do.  There may already be a bug reported about it.

Octave could also switch to using the gnulib module for log2.  This would also fix the error, but only on a native build.

Rik <rik5>
Group administrator
Thu 19 Jun 2014 05:17:02 PM UTC, comment #6: 

So it occurs that the problem is actually in the Cygwin C library which implements log2 in terms of log (I discover at this occasion that log2 has been standardized in C99 and C++11). Should I send a bug report to the Cygwin mailing list?

Falk Tannhäuser <fata>
Thu 19 Jun 2014 03:24:26 PM UTC, comment #5: 

This is likely to be a problem with 1) cross-compiling, 2) the availibility of a log2 function in the C library, or 3) Cygwin's implementation of log2.

For efficiency, Octave tries to map log2 to the local C library function which is usually a better choice than trying to calculate it ourselves.  The code is in liboctave/numeric/lo-mappers.cc.


double
xlog2 (double x)
{
#if defined (HAVE_LOG2)
  return log2 (x);
#else
#if defined (M_LN2)
  static double ln2 = M_LN2;
#else
  static double ln2 = log (2);
#endif

  return log (x) / ln2;
#endif
}


If I check my config.h file I see that HAVE_LOG2 is defined to 1 so glibc has the log2 function.

How did you get Octave for cygwin?  Is this a native build or a cross-compiled build?  If it is a cross-compiled build using the MXE framework then it may be that the tools could not determine whether the end system would have log2 in libc and so undefined HAVE_LOG2.

Another possibility is that the log2 function on cygwin is just broken.  The gnulib project has replacement functions for ones that are known to have compatibility problems between platforms.  If I look at log2 (https://www.gnu.org/software/gnulib/manual/html_node/log2.html) I see that Cygwin is listed.


This function returns slightly wrong values for exact powers of 2 on some platforms: Cygwin 1.7.9.


Rik <rik5>
Group administrator
Thu 19 Jun 2014 08:54:39 AM UTC, comment #4: 

With Octave 3.8.2 compiled with gcc 4.7 0n OSX 10.9.3,
I get:


>> n = -1074:1023;
>> errors = find(log2(pow2(n)) ~= n) - 1075
errors = []



Carlo de Falco <cdf>
Group Member
Thu 19 Jun 2014 06:30:07 AM UTC, comment #3: 

Interestingly, with

n = -323:308;
errors = find(log10(10.^n) ~= n) - 324


I get

errors = -323 : -312

i.e. despite the fact that 10.^n can't be precisely represented for n<0 or n>22, errors occur only for denormalized arguments, which is unsurprising.

Falk Tannhäuser <fata>
Thu 19 Jun 2014 06:05:13 AM UTC, comment #2: 

I have Octave 3.8.1 configured for "x86_64-unknown-cygwin".
octave_config_info.GCC_VERSION is 4.8.2.

Falk Tannhäuser <fata>
Thu 19 Jun 2014 01:53:06 AM UTC, comment #1: 

This works fine for me, i.e., errors is empty on a Linux build with gcc 4.6.3.

Which version of Octave for Windows are you using?

Rik <rik5>
Group administrator
Wed 18 Jun 2014 11:59:41 PM UTC, original submission:  

2^n is exactly representable as double for integer n between -1074 and 1023. However, the following code


n = -1074:1023;
errors = find(log2(pow2(n)) ~= n) - 1075


shows that log2 returns inaccurate results for 441 values of n (among others: -1066, -1023,  -1021, ..., -39, -31, -29, 29, 31, 39, ..., 1019, 1021, 1023).
OTOH, with Matlab version 8.0.0.783 (R2012b) an accurate result is always returned and the "errors" variable is empty with above code.

By the way, the FPU of Intel x86 has an instruction "FYL2X" giving accurate results if 1.0 is loaded for Y argument, but unfortunately the C/C++ standard math library has only log() and log10() functions. Rounding errors occur if log2() is implemented in terms of log(), as can be seen with the following C++ program:


#include <cmath>
#include <limits>
#include <iostream>

#if 1 //////////////////////////////////////////////////////////
inline double log2(double x) { return std::log(x)/std::log(2); }
#else //////////////////////////////////////////////////////////
inline double log2(double x)
{
  double one = 1;
  double result;
  asm ("fyl2x" : "=t" (result) : "0" (x), "u" (one) : "st(1)");
  return result;
}
#endif /////////////////////////////////////////////////////////

int main()
{
  int n = 0;
  for(double x = 1; x < std::numeric_limits<double>::infinity(); x*=2, ++n)
  //for(double x = 1; x > 0; x/=2, --n)
    if(n != log2(x))
      std::cout << n << '\n';
  return 0;
}


The assembler version of log2() compiles e.g. with GCC 4.9.0 for target x86_64-pc-cygwin.

Falk Tannhäuser <fata>

 

(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 dasergatskov (Posted a comment)
  • -email is unavailable- added by dac922 (Posted a comment)
  • -email is unavailable- added by siko1056 (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by cdf (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by fata (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 4 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2014-06-22 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2014-06-19 rik5 StatusWorks For Me Ready For Test
    2014-06-19 rik5 StatusNone Works For Me

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code