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

 
 

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

Submitted by:  Falk Tannhäuser <fata>
Submitted on:  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

Add a New Comment (Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

( Jump to the original submission)

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>
Project 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>
Project 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>
Project 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>
Project 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>
Project 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>
Project 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>
Project 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>
Project 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>
Project 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>
Project 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>
Project 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>
Project 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 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)
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project 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.7