octave-bug-tracker
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Octave-bug-tracker] [bug #42583] log2() returns inaccurate result for m


From: Falk Tannhäuser
Subject: [Octave-bug-tracker] [bug #42583] log2() returns inaccurate result for many integer powers of 2, unlike Matlab
Date: Wed, 18 Jun 2014 23:59:41 +0000
User-agent: Mozilla/5.0 (Windows NT 6.3; WOW64; rv:30.0) Gecko/20100101 Firefox/30.0

URL:
  <http://savannah.gnu.org/bugs/?42583>

                 Summary: log2() returns inaccurate result for many integer
powers of 2, unlike Matlab
                 Project: GNU Octave
            Submitted by: fata
            Submitted on: Do 19 Jun 2014 01:59:41 CEST
                Category: Octave Function
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Inaccurate Result
                  Status: None
             Assigned to: None
         Originator Name: Falk Tannhäuser
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any
                 Release: 3.8.1
        Operating System: Microsoft Windows

    _______________________________________________________

Details:

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.




    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?42583>

_______________________________________________
  Nachricht gesendet von/durch Savannah
  http://savannah.gnu.org/




reply via email to

[Prev in Thread] Current Thread [Next in Thread]