log2() function from C standard library is inaccurate

[Marco Atzeri]

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