diff options
| author | Jeff Johnston <jjohnstn@redhat.com> | 2010-03-08 20:16:37 +0300 |
|---|---|---|
| committer | Jeff Johnston <jjohnstn@redhat.com> | 2010-03-08 20:16:37 +0300 |
| commit | 23a6adc2c38454ff01231a82841e9c8a0c09ff45 (patch) | |
| tree | da7cb975677a93aca19967efbd3308dd623b3e54 | |
| parent | 1f440e4dbe41b1237c7a75c0723506a5cc2003d3 (diff) | |
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.
| -rw-r--r-- | newlib/ChangeLog | 9 | ||||
| -rw-r--r-- | newlib/libm/common/s_rint.c | 30 | ||||
| -rw-r--r-- | newlib/libm/common/sf_rint.c | 2 |
3 files changed, 29 insertions, 12 deletions
diff --git a/newlib/ChangeLog b/newlib/ChangeLog index e2791eacf..578e297d9 100644 --- a/newlib/ChangeLog +++ b/newlib/ChangeLog @@ -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. diff --git a/newlib/libm/common/s_rint.c b/newlib/libm/common/s_rint.c index 4fa5ebc6c..76cff08a1 100644 --- a/newlib/libm/common/s_rint.c +++ b/newlib/libm/common/s_rint.c @@ -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 */ - - - diff --git a/newlib/libm/common/sf_rint.c b/newlib/libm/common/sf_rint.c index 6459b7a4c..bc0b46659 100644 --- a/newlib/libm/common/sf_rint.c +++ b/newlib/libm/common/sf_rint.c @@ -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 */ |