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@dessent.net>
* 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;