mirror of
git://sourceware.org/git/newlib-cygwin.git
synced 2025-01-23 23:47:22 +08:00
9e329b544a
soft-fp should round floating pointer numbers according to the current
rounding mode. However, in the current code of lrint() and llrint(),
there are if statements before the actual rounding computation
if(j0 < -1)
return 0;
Where j0 is the exponent of the floating point number.
It means any number having a exponent less than -1
(i.e. interval (-0.5, 0.5)) will be rounded to 0 regardeless of the
rounding mode.
The bug already fixed in glibc in 2006 by moving the check afterwards
the rounding computation, but still persists in newlib.
This patch fixed it in a similar way to glibc
Ref Commit in glibc: 6624dbc07b5a9fb316ed188ef01f65b8eea8b47c
116 lines
3.0 KiB
C
116 lines
3.0 KiB
C
/* lrint adapted to be llrint for Newlib, 2009 by Craig Howland. */
|
|
/* @(#)s_lrint.c 5.1 93/09/24 */
|
|
/*
|
|
* ====================================================
|
|
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
|
*
|
|
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
|
* Permission to use, copy, modify, and distribute this
|
|
* software is freely granted, provided that this notice
|
|
* is preserved.
|
|
* ====================================================
|
|
*/
|
|
|
|
/*
|
|
* llrint(x)
|
|
* Return x rounded to integral value according to the prevailing
|
|
* rounding mode.
|
|
* Method:
|
|
* Using floating addition.
|
|
* Exception:
|
|
* Inexact flag raised if x not equal to llrint(x).
|
|
*/
|
|
|
|
#include "fdlibm.h"
|
|
|
|
#ifndef _DOUBLE_IS_32BITS
|
|
|
|
#ifdef __STDC__
|
|
static const double
|
|
#else
|
|
static double
|
|
#endif
|
|
|
|
/* Adding a double, x, to 2^52 will cause the result to be rounded based on
|
|
the fractional part of x, according to the implementation's current rounding
|
|
mode. 2^52 is the smallest double that can be represented using all 52 significant
|
|
digits. */
|
|
TWO52[2]={
|
|
4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
|
|
-4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
|
|
};
|
|
|
|
long long int
|
|
#ifdef __STDC__
|
|
llrint(double x)
|
|
#else
|
|
llrint(x)
|
|
double x;
|
|
#endif
|
|
{
|
|
__int32_t i0,j0,sx;
|
|
__uint32_t i1;
|
|
double t;
|
|
volatile double w;
|
|
long long int result;
|
|
|
|
EXTRACT_WORDS(i0,i1,x);
|
|
|
|
/* Extract sign bit. */
|
|
sx = (i0>>31)&1;
|
|
|
|
/* Extract exponent field. */
|
|
j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
|
|
/* j0 in [-1023,1024] */
|
|
|
|
if(j0 < 20)
|
|
{
|
|
/* j0 in [-1023,19] */
|
|
w = TWO52[sx] + x;
|
|
t = w - TWO52[sx];
|
|
GET_HIGH_WORD(i0, t);
|
|
/* Detect the all-zeros representation of plus and
|
|
minus zero, which fails the calculation below. */
|
|
if ((i0 & ~((__int32_t)1 << 31)) == 0)
|
|
return 0;
|
|
j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
|
|
i0 &= 0x000fffff;
|
|
i0 |= 0x00100000;
|
|
result = (j0 < 0 ? 0 : i0 >> (20 - j0));
|
|
}
|
|
else if (j0 < (int)(8 * sizeof (long long int)) - 1)
|
|
{
|
|
/* 64bit return: j0 in [20,62] */
|
|
if (j0 >= 52)
|
|
/* 64bit return: j0 in [52,62] */
|
|
/* 64bit return: left shift amt in [32,42] */
|
|
result = ((long long int) ((i0 & 0x000fffff) | 0x00100000) << (j0 - 20)) |
|
|
/* 64bit return: right shift amt in [0,10] */
|
|
((long long int) i1 << (j0 - 52));
|
|
else
|
|
{
|
|
/* 64bit return: j0 in [20,51] */
|
|
w = TWO52[sx] + x;
|
|
t = w - TWO52[sx];
|
|
EXTRACT_WORDS (i0, i1, t);
|
|
j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
|
|
i0 &= 0x000fffff;
|
|
i0 |= 0x00100000;
|
|
/* After round:
|
|
* 64bit return: j0 in [20,52] */
|
|
/* 64bit return: left shift amt in [0,32] */
|
|
/* ***64bit return: right shift amt in [32,0] */
|
|
result = ((long long int) i0 << (j0 - 20))
|
|
| SAFE_RIGHT_SHIFT (i1, (52 - j0));
|
|
}
|
|
}
|
|
else
|
|
{
|
|
return (long long int) x;
|
|
}
|
|
|
|
return sx ? -result : result;
|
|
}
|
|
|
|
#endif /* _DOUBLE_IS_32BITS */
|