This is the mail archive of the
cygwin
mailing list for the Cygwin project.
Re: log2() function from C standard library is inaccurate
- From: Marco Atzeri <marco dot atzeri at gmail dot com>
- To: cygwin at cygwin dot com
- Date: Fri, 20 Jun 2014 09:55:39 +0200
- Subject: Re: log2() function from C standard library is inaccurate
- Authentication-results: sourceware.org; auth=none
- References: <53A3715B dot 8000904 at free dot fr>
On 20/06/2014 01:25, Falk Tannhäuser wrote:
I noticed that the log2() function in GNU Octave returns inaccurate results
for many arguments that are exactly representable integer powers of 2.
After discussion on the Octave mailing list
https://savannah.gnu.org/bugs/?42583
it occurs that Octave doesn't have this problem on most other platforms
(Linux, OS X) and that it is due to the fact that Octave usually uses
log2() from the local C library.
The following C++ program exposes the error:
_________________________________________________________________
// Compilation: g++ -Wall -Wextra -O3 -s testlog2.cxx -o testlog2
#include <limits>
#include <iostream>
#if 0 /////////////////////////////////////////////////////////
inline double log2(double x)
{
double result;
asm ("fyl2x" : "=t" (result) : "0" (x), "u" (1.0) : "st(1)");
return result;
}
#else /////////////////////////////////////////////////////////
#include <math.h>
#endif ////////////////////////////////////////////////////////
int main()
{
int n = std::numeric_limits<double>::min_exponent
- std::numeric_limits<double>::digits;
std::cout << "n = " << n << " ... "
<< std::numeric_limits<double>::max_exponent - 1
<< '\n';
int errcnt = 0;
for(double x = std::numeric_limits<double>::denorm_min();
x < std::numeric_limits<double>::infinity();
x *= 2, ++n)
{
if(n != log2(x))
{
std::cout << n << '\n';
++errcnt;
}
}
std::cout << "Found " << errcnt << " errors.\n";
return 0;
}
_________________________________________________________________
With GCC 4.9.0 or 4.8.3 for x86_64-pc-cygwin, erroneous results
are found for 441 values.
for 32 bit 2075, due to the incorrect comparison.
This seems to be due to the fact that
log2 is implemented in terms of log (as something like log(x)/log(2)).
OTOH, when using the assembler version of log2 (just changing '#if 0'
to '#if 1' in above code); accurate results are always returned.
Your analysis is correct,
from the source "newlib/libm/common/s_log2.c"
-----------------------------------------------------
The Newlib implementations are not full, intrinisic calculations, but
rather are derivatives based on <<log>>. (Accuracy might be slightly
off from
a direct calculation.) In addition to functions, they are also
implemented as
macros defined in math.h:
. #define log2(x) (log (x) / _M_LN2)
. #define log2f(x) (logf (x) / (float) _M_LN2)
------------------------------------------------------
Falk
The C library project used by cygwin is here
https://sourceware.org/newlib/
Regards
Marco
--
Problem reports: http://cygwin.com/problems.html
FAQ: http://cygwin.com/faq/
Documentation: http://cygwin.com/docs.html
Unsubscribe info: http://cygwin.com/ml/#unsubscribe-simple