X-Recipient: archive-cygwin AT delorie DOT com
DomainKey-Signature: a=rsa-sha1; c=nofws; d=sourceware.org; h=list-id
	:list-unsubscribe:list-subscribe:list-archive:list-post
	:list-help:sender:message-id:date:from:mime-version:to:subject
	:references:in-reply-to:content-type:content-transfer-encoding;
	 q=dns; s=default; b=tvmuGwBB4cPG0piCPAzcqRiZaf3fRQ1da8Rk6HRA+9z
	HZuB9PFcdWKZTCzwLKFAHtOyT4SYwhkwFmilNq5nb8dnlCcbEj7wUKe7F2CHW1pd
	3lUi6Z8OGlpZKuqWXiU1EzT7wnQJ4Rtn1xRexOQRWbX7/y6GcD/s+OAI4/VvrxAA
	=
DKIM-Signature: v=1; a=rsa-sha1; c=relaxed; d=sourceware.org; h=list-id
	:list-unsubscribe:list-subscribe:list-archive:list-post
	:list-help:sender:message-id:date:from:mime-version:to:subject
	:references:in-reply-to:content-type:content-transfer-encoding;
	 s=default; bh=mrGAADc6EZh9mnjZDQAAd2e37WQ=; b=DUSalcI8ppHkYCaw4
	kUY+PEwOHjG1tEDtrHYYkEGSdyEpGrlxkqNZZS85alcJ8yp9i7PB5Z2hWXpch60w
	QganP6yUPHpTk0VbOKBtTMSugV/3HCFwJQg66sTUzpBiCSajKSoUQ61pwxf6kxV9
	wpC9eX3r+peAJxm1ndpKxX3IPg=
Mailing-List: contact cygwin-help AT cygwin DOT com; run by ezmlm
List-Id: <cygwin.cygwin.com>
List-Subscribe: <mailto:cygwin-subscribe AT cygwin DOT com>
List-Archive: <http://sourceware.org/ml/cygwin/>
List-Post: <mailto:cygwin AT cygwin DOT com>
List-Help: <mailto:cygwin-help AT cygwin DOT com>, <http://sourceware.org/ml/#faqs>
Sender: cygwin-owner AT cygwin DOT com
Mail-Followup-To: cygwin AT cygwin DOT com
Delivered-To: mailing list cygwin AT cygwin DOT com
Authentication-Results: sourceware.org; auth=none
X-Virus-Found: No
X-Spam-SWARE-Status: No, score=-1.9 required=5.0 tests=AWL,BAYES_00,FREEMAIL_FROM,RCVD_IN_DNSWL_LOW,SPF_PASS autolearn=ham version=3.3.2
X-HELO: mail-wi0-f178.google.com
X-Received: by 10.180.94.37 with SMTP id cz5mr2122776wib.19.1403250942879;        Fri, 20 Jun 2014 00:55:42 -0700 (PDT)
Message-ID: <53A3E8FB.7030702@gmail.com>
Date: Fri, 20 Jun 2014 09:55:39 +0200
From: Marco Atzeri <marco DOT atzeri AT gmail DOT com>
User-Agent: Mozilla/5.0 (Windows NT 6.1; WOW64; rv:24.0) Gecko/20100101 Thunderbird/24.6.0
MIME-Version: 1.0
To: cygwin AT cygwin DOT com
Subject: Re: log2() function from C standard library is inaccurate
References: <53A3715B DOT 8000904 AT free DOT fr>
In-Reply-To: <53A3715B.8000904@free.fr>
Content-Type: text/plain; charset=ISO-8859-15; format=flowed
Content-Transfer-Encoding: 8bit
X-IsSubscribed: yes

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