JehanneOS
/
newlib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

2010-03-08 Craig Howland <howland@LGSInnovations.com> 

* libm/common/s_rint.c:  Fix error when integral part had 18 bits and
        fraction had bits set beyond first radix bit.  Also, make 2-part
        adjustment consistent with 1-part adjustment when adjusting fractional
        bits.
        * libm/common/sf_rint.c:  Make fractional-bit adjustment consistent
        with s_rint.c by setting 0b.01 instead of 0b.001.
This commit is contained in:
Jeff Johnston 2010-03-08 17:16:37 +00:00
parent 1f440e4dbe
commit 23a6adc2c3
3 changed files with 29 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
newlib/ChangeLog
View File

@ -1,3 +1,12 @@
2010-03-08 Craig Howland <howland@LGSInnovations.com>
* libm/common/s_rint.c: Fix error when integral part had 18 bits and
fraction had bits set beyond first radix bit. Also, make 2-part
adjustment consistent with 1-part adjustment when adjusting fractional
bits.
* libm/common/sf_rint.c: Make fractional-bit adjustment consistent
with s_rint.c by setting 0b.01 instead of 0b.001.
2010-03-05 Craig Howland <howland@LGSInnovations.com>
* libm/math/ef_sqrt.c: Delete unused variable sign.

30
newlib/libm/common/s_rint.c
View File

@ -51,6 +51,16 @@ SEEALSO
* rounding mode.
* Method:
* Using floating addition.
* Whenever a fraction is present, if the second or any following bit after
* the radix point is set, limit to the second radix point to avoid
* possible double rounding in the TWO52 +- steps (in case guard bits are
* used). Specifically, if have any, chop off bits past the 2nd place and
* set the second place.
* (e.g. 2.0625=0b10.0001 => 0b10.01=2.25;
* 2.3125=0b10.011 => 0b10.01=2.25;
* 1.5625= 0b1.1001 => 0b1.11=1.75;
* 1.9375= 0b1.1111 => 0b1.11=1.75.
* Pseudo-code: if(x.frac & ~0b0.10) x.frac = (x.frac & 0b0.11) | 0b0.01;).
* Exception:
* Inexact flag raised if x not equal to rint(x).
*/
@ -81,11 +91,11 @@ TWO52[2]={
double t;
volatile double w;
EXTRACT_WORDS(i0,i1,x);
sx = (i0>>31)&1;
j0 = ((i0>>20)&0x7ff)-0x3ff;
if(j0<20) {
if(j0<0) {
if(((i0&0x7fffffff)|i1)==0) return x;
sx = (i0>>31)&1; /* sign */
j0 = ((i0>>20)&0x7ff)-0x3ff; /* exponent */
if(j0<20) { /* no integral bits in LS part */
if(j0<0) { /* x is fractional or 0 */
if(((i0&0x7fffffff)|i1)==0) return x; /* x == 0 */
i1 |= (i0&0x0fffff);
i0 &= 0xfffe0000;
i0 |= ((i1|-i1)>>12)&0x80000;
@ -95,13 +105,14 @@ TWO52[2]={
GET_HIGH_WORD(i0,t);
SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
return t;
} else {
} else { /* x has integer and maybe fraction */
i = (0x000fffff)>>j0;
if(((i0&i)|i1)==0) return x; /* x is integral */
i>>=1;
if(((i0&i)|i1)!=0) {
if(j0==19) i1 = 0x40000000; else
i0 = (i0&(~i))|((0x20000)>>j0);
/* 2nd or any later bit after radix is set */
if(j0==19) i1 = 0x80000000; else i1 = 0;
i0 = (i0&(~i))|((0x40000)>>j0);
}
}
} else if (j0>51) {
@ -119,6 +130,3 @@ TWO52[2]={
}
#endif /* _DOUBLE_IS_32BITS */

2
newlib/libm/common/sf_rint.c
View File

@ -57,7 +57,7 @@ TWO23[2]={
i = (0x007fffff)>>j0;
if((i0&i)==0) return x; /* x is integral */
i>>=1;
if((i0&i)!=0) i0 = (i0&(~i))|((0x100000)>>j0);
if((i0&i)!=0) i0 = (i0&(~i))|((0x200000)>>j0);
}
} else {
if(!FLT_UWORD_IS_FINITE(ix)) return x+x; /* inf or NaN */