diff options
Diffstat (limited to 'newlib/libm/common/s_rint.c')
| -rw-r--r-- | newlib/libm/common/s_rint.c | 30 |
1 files changed, 19 insertions, 11 deletions
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 */ - - - |