Subject: Re: [patch] lrint/lrintf oddity
From: Jeff Johnston <[email-masked]>
To: Brian Dessent <[email-masked]>
Cc: [email-masked], [email-masked]
Date: Thu, 08 Sep 2005 17:05:19 -0400

Thanks Brian.  Patch checked in and s_lrint.c modified similiarly.

-- Jeff J.

Brian Dessent wrote:
> 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  <[email-masked]>
> 
> 	* sf_lrint.c (lrintf): Mask 'i0' correctly when extracting
> 	mantissa.
> 
> 
> ------------------------------------------------------------------------
> 
> 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;
> 
> 

--
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/