Mailing-List: contact cygwin-help AT cygwin DOT com; run by ezmlm
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
Message-ID: <4319C016.2958BC77@dessent.net>
Date: Sat, 03 Sep 2005 08:24:06 -0700
From: Brian Dessent <brian AT dessent DOT net>
MIME-Version: 1.0
To: cygwin AT cygwin DOT com, newlib AT sources DOT redhat DOT com
Subject: [patch] lrint/lrintf oddity
References: <20050903185145 DOT 5bb6c57e DOT cygwin-erikd AT mega-nerd DOT com>
Content-Type: multipart/mixed;  boundary="------------7CF6A2CE1526A4D41737FB5D"
X-IsSubscribed: yes
Reply-To: cygwin AT cygwin DOT com

--------------7CF6A2CE1526A4D41737FB5D
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

Erik de Castro Lopo wrote:

>     $ gcc -Wall -O3 lrintf_test.c -o lrintf_test
>     $ ./lrintf_test.exe
>     lrintf_test ( 3333.000000) : ok
>     lrint_test  ( 3333.000000) : ok
>     lrintf_test (33333.000000) :
> 
>     Line 70: float : Incorrect sample (#252 : -8399916.000000 => 889181140).

Yes, this does in fact seem to be a bug in lrintf.  A reduced testcase
is:

#include <stdio.h>
#include <math.h>

int main ()
{
  float f = -8399916.0;
  printf ("f = %f, k = %ld\n", f, lrintf (f));
}

The value -8399916.0 in binary is

1 10010110 00000000010110000101100

The exponent (j0) is 23 (150 - 127 = 23).

In lrintf, the problem is the part of the algorithm that acts on the
case of 23 <= j0 < 31.  In this case, the mantissa is simply extracted
and shifted left by (j0 - 23).  This is handled by line 66 of
sf_lrint.c:

        result = (long int) i0 << (j0 - 23);

This variable 'result' is then returned, optionally negated if the sign
bit was set.  However, this calculation makes two errors: First, it must
mask off only the mantissa from i0, removing the sign bit and exponent. 
Second, it neglects to add the implicit MSB "1" which is always present
to the left of the decimal but not actually included in the mantissa.

The attached patch fixes these two problems and fixes the poster's test
case. It appears that lrint() also might suffer from the same problem:

  else if (j0 < (int)(8 * sizeof (long int)) - 1)
    {
      if (j0 >= 52)
        result = ((long int) i0 << (j0 - 20)) | (i1 << (j0 - 52));

However, this section of code will never be reached on 32 bit hardware
since the first 'if' will guarantee that j0 < 31.  However this looks
like it's still broken for 64 bit machines, and someone should probably
audit this.

Brian

2005-09-03  Brian Dessent  <brian AT dessent DOT net>

	* sf_lrint.c (lrintf): Mask 'i0' correctly when extracting
	mantissa.
--------------7CF6A2CE1526A4D41737FB5D
Content-Type: text/plain; charset=us-ascii;
 name="newlib-lrintf.patch"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="newlib-lrintf.patch"

Index: sf_lrint.c
===================================================================
RCS file: /cvs/src/src/newlib/libm/common/sf_lrint.c,v
retrieving revision 1.2
diff -u -p -r1.2 sf_lrint.c
--- sf_lrint.c	28 Jun 2005 17:03:18 -0000	1.2
+++ sf_lrint.c	3 Sep 2005 15:13:06 -0000
@@ -63,7 +63,7 @@ TWO23[2]={
       if (j0 < -1)
         return 0;
       else if (j0 >= 23)
-        result = (long int) i0 << (j0 - 23);
+        result = (long int) (i0 & 0x7fffff | 0x800000) << (j0 - 23);
       else
         {
           w = TWO23[sx] + x;




--------------7CF6A2CE1526A4D41737FB5D
Content-Type: text/plain; charset=us-ascii

--
Unsubscribe info:      http://cygwin.com/ml/#unsubscribe-simple
Problem reports:       http://cygwin.com/problems.html
Documentation:         http://cygwin.com/docs.html
FAQ:                   http://cygwin.com/faq/
--------------7CF6A2CE1526A4D41737FB5D--