1
0
mirror of https://git.FreeBSD.org/src.git synced 2025-01-20 15:43:16 +00:00

Fix a bug where rintf() rounded the wrong way in round-to-nearest mode

on all inputs of the form x.75, where x is an even integer and
log2(x) = 21.  A similar problem occurred when rounding upward.
The bug involves the following snippet copied from rint():

	i>>=1;
	if((i0&i)!=0) i0 = (i0&(~i))|((0x100000)>>j0);

The constant 0x100000 should be 0x200000.  Apparently this case was
never tested.

It turns out that the bit manipulation is completely superfluous
anyway, so remove it.  (It tries to simulate 90% of the rounding
process that the FPU does anyway.)  Also, the special case of +-0 is
handled twice (in different ways), so remove the second instance.

Throw in some related simplifications from bde:

- Work around a bug where gcc fails to clip to float precision by
  declaring two float variables as volatile.  Previously, we
  tricked gcc into generating correct code by declaring some
  float constants as doubles.

- Remove additional superfluous bit manipulation.

- Minor reorganization.

- Include <sys/types.h> explicitly.

Note that some of the equivalent lines in rint() also appear to be
unnecessary, but I'll defer to the numerical analysts who wrote it,
since I can't test all 2^64 cases.

Discussed with:	bde
This commit is contained in:
David Schultz 2004-06-09 21:24:52 +00:00
parent 782de70aff
commit c4da2324a3
Notes: svn2git 2020-12-20 02:59:44 +00:00
svn path=/head/; revision=130285

View File

@ -17,16 +17,11 @@
static char rcsid[] = "$FreeBSD$";
#endif
#include <sys/types.h>
#include "math.h"
#include "math_private.h"
/*
* TWO23 is double instead of float to avoid a bug in gcc. Without
* this, gcc thinks that TWO23[sx]+x and w-TWO23[sx] already have float
* precision and doesn't clip them to float precision when they are
* assigned and returned.
*/
static const double
static const float
TWO23[2]={
8.3886080000e+06, /* 0x4b000000 */
-8.3886080000e+06, /* 0xcb000000 */
@ -36,34 +31,20 @@ float
rintf(float x)
{
int32_t i0,j0,sx;
u_int32_t i,i1;
float w,t;
volatile float w,t; /* volatile works around gcc bug */
GET_FLOAT_WORD(i0,x);
sx = (i0>>31)&1;
j0 = ((i0>>23)&0xff)-0x7f;
if(j0<23) {
if(j0<0) {
if((i0&0x7fffffff)==0) return x;
i1 = (i0&0x07fffff);
i0 &= 0xfff00000;
i0 |= ((i1|-i1)>>9)&0x400000;
SET_FLOAT_WORD(x,i0);
w = TWO23[sx]+x;
t = w-TWO23[sx];
GET_FLOAT_WORD(i0,t);
SET_FLOAT_WORD(t,(i0&0x7fffffff)|(sx<<31));
return t;
} else {
i = (0x007fffff)>>j0;
if((i0&i)==0) return x; /* x is integral */
i>>=1;
if((i0&i)!=0) i0 = (i0&(~i))|((0x100000)>>j0);
}
} else {
if(j0==0x80) return x+x; /* inf or NaN */
else return x; /* x is integral */
w = TWO23[sx]+x;
return w-TWO23[sx];
}
SET_FLOAT_WORD(x,i0);
w = TWO23[sx]+x;
return w-TWO23[sx];
if(j0==0x80) return x+x; /* inf or NaN */
else return x; /* x is integral */
}