diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/e_remainder.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_remainder.c | 180 |
1 files changed, 101 insertions, 79 deletions
diff --git a/sysdeps/ieee754/dbl-64/e_remainder.c b/sysdeps/ieee754/dbl-64/e_remainder.c index 2d20bb1dfe..c6a1901b10 100644 --- a/sysdeps/ieee754/dbl-64/e_remainder.c +++ b/sysdeps/ieee754/dbl-64/e_remainder.c @@ -39,89 +39,111 @@ /* An ultimate remainder routine. Given two IEEE double machine numbers x */ /* ,y it computes the correctly rounded (to nearest) value of remainder */ /**************************************************************************/ -double __ieee754_remainder(double x, double y) +double +__ieee754_remainder (double x, double y) { - double z,d,xx; - int4 kx,ky,n,nn,n1,m1,l; - mynumber u,t,w={{0,0}},v={{0,0}},ww={{0,0}},r; - u.x=x; - t.x=y; - kx=u.i[HIGH_HALF]&0x7fffffff; /* no sign for x*/ - t.i[HIGH_HALF]&=0x7fffffff; /*no sign for y */ - ky=t.i[HIGH_HALF]; + double z, d, xx; + int4 kx, ky, n, nn, n1, m1, l; + mynumber u, t, w = { { 0, 0 } }, v = { { 0, 0 } }, ww = { { 0, 0 } }, r; + u.x = x; + t.x = y; + kx = u.i[HIGH_HALF] & 0x7fffffff; /* no sign for x*/ + t.i[HIGH_HALF] &= 0x7fffffff; /*no sign for y */ + ky = t.i[HIGH_HALF]; /*------ |x| < 2^1023 and 2^-970 < |y| < 2^1024 ------------------*/ - if (kx<0x7fe00000 && ky<0x7ff00000 && ky>=0x03500000) { - SET_RESTORE_ROUND_NOEX (FE_TONEAREST); - if (kx+0x00100000<ky) return x; - if ((kx-0x01500000)<ky) { - z=x/t.x; - v.i[HIGH_HALF]=t.i[HIGH_HALF]; - d=(z+big.x)-big.x; - xx=(x-d*v.x)-d*(t.x-v.x); - if (d-z!=0.5&&d-z!=-0.5) return (xx!=0)?xx:((x>0)?ZERO.x:nZERO.x); - else { - if (ABS(xx)>0.5*t.x) return (z>d)?xx-t.x:xx+t.x; - else return xx; - } - } /* (kx<(ky+0x01500000)) */ - else { - r.x=1.0/t.x; - n=t.i[HIGH_HALF]; - nn=(n&0x7ff00000)+0x01400000; - w.i[HIGH_HALF]=n; - ww.x=t.x-w.x; - l=(kx-nn)&0xfff00000; - n1=ww.i[HIGH_HALF]; - m1=r.i[HIGH_HALF]; - while (l>0) { - r.i[HIGH_HALF]=m1-l; - z=u.x*r.x; - w.i[HIGH_HALF]=n+l; - ww.i[HIGH_HALF]=(n1)?n1+l:n1; - d=(z+big.x)-big.x; - u.x=(u.x-d*w.x)-d*ww.x; - l=(u.i[HIGH_HALF]&0x7ff00000)-nn; - } - r.i[HIGH_HALF]=m1; - w.i[HIGH_HALF]=n; - ww.i[HIGH_HALF]=n1; - z=u.x*r.x; - d=(z+big.x)-big.x; - u.x=(u.x-d*w.x)-d*ww.x; - if (ABS(u.x)<0.5*t.x) return (u.x!=0)?u.x:((x>0)?ZERO.x:nZERO.x); + if (kx < 0x7fe00000 && ky < 0x7ff00000 && ky >= 0x03500000) + { + SET_RESTORE_ROUND_NOEX (FE_TONEAREST); + if (kx + 0x00100000 < ky) + return x; + if ((kx - 0x01500000) < ky) + { + z = x / t.x; + v.i[HIGH_HALF] = t.i[HIGH_HALF]; + d = (z + big.x) - big.x; + xx = (x - d * v.x) - d * (t.x - v.x); + if (d - z != 0.5 && d - z != -0.5) + return (xx != 0) ? xx : ((x > 0) ? ZERO.x : nZERO.x); + else + { + if (ABS (xx) > 0.5 * t.x) + return (z > d) ? xx - t.x : xx + t.x; + else + return xx; + } + } /* (kx<(ky+0x01500000)) */ else - if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x; - else - {z=u.x/t.x; d=(z+big.x)-big.x; return ((u.x-d*w.x)-d*ww.x);} - } - - } /* (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000) */ - else { - if (kx<0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) { - y=ABS(y)*t128.x; - z=__ieee754_remainder(x,y)*t128.x; - z=__ieee754_remainder(z,y)*tm128.x; - return z; - } - else { - if ((kx&0x7ff00000)==0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) { - y=ABS(y); - z=2.0*__ieee754_remainder(0.5*x,y); - d = ABS(z); - if (d <= ABS(d-y)) return z; - else return (z>0)?z-y:z+y; - } - else { /* if x is too big */ - if (ky==0 && t.i[LOW_HALF] == 0) /* y = 0 */ - return (x * y) / (x * y); - else if (kx >= 0x7ff00000 /* x not finite */ - || (ky>0x7ff00000 /* y is NaN */ - || (ky == 0x7ff00000 && t.i[LOW_HALF] != 0))) - return (x * y) / (x * y); + { + r.x = 1.0 / t.x; + n = t.i[HIGH_HALF]; + nn = (n & 0x7ff00000) + 0x01400000; + w.i[HIGH_HALF] = n; + ww.x = t.x - w.x; + l = (kx - nn) & 0xfff00000; + n1 = ww.i[HIGH_HALF]; + m1 = r.i[HIGH_HALF]; + while (l > 0) + { + r.i[HIGH_HALF] = m1 - l; + z = u.x * r.x; + w.i[HIGH_HALF] = n + l; + ww.i[HIGH_HALF] = (n1) ? n1 + l : n1; + d = (z + big.x) - big.x; + u.x = (u.x - d * w.x) - d * ww.x; + l = (u.i[HIGH_HALF] & 0x7ff00000) - nn; + } + r.i[HIGH_HALF] = m1; + w.i[HIGH_HALF] = n; + ww.i[HIGH_HALF] = n1; + z = u.x * r.x; + d = (z + big.x) - big.x; + u.x = (u.x - d * w.x) - d * ww.x; + if (ABS (u.x) < 0.5 * t.x) + return (u.x != 0) ? u.x : ((x > 0) ? ZERO.x : nZERO.x); + else + if (ABS (u.x) > 0.5 * t.x) + return (d > z) ? u.x + t.x : u.x - t.x; + else + { + z = u.x / t.x; d = (z + big.x) - big.x; + return ((u.x - d * w.x) - d * ww.x); + } + } + } /* (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000) */ + else + { + if (kx < 0x7fe00000 && ky < 0x7ff00000 && (ky > 0 || t.i[LOW_HALF] != 0)) + { + y = ABS (y) * t128.x; + z = __ieee754_remainder (x, y) * t128.x; + z = __ieee754_remainder (z, y) * tm128.x; + return z; + } else - return x; + { + if ((kx & 0x7ff00000) == 0x7fe00000 && ky < 0x7ff00000 && + (ky > 0 || t.i[LOW_HALF] != 0)) + { + y = ABS (y); + z = 2.0 * __ieee754_remainder (0.5 * x, y); + d = ABS (z); + if (d <= ABS (d - y)) + return z; + else + return (z > 0) ? z - y : z + y; + } + else /* if x is too big */ + { + if (ky == 0 && t.i[LOW_HALF] == 0) /* y = 0 */ + return (x * y) / (x * y); + else if (kx >= 0x7ff00000 /* x not finite */ + || (ky > 0x7ff00000 /* y is NaN */ + || (ky == 0x7ff00000 && t.i[LOW_HALF] != 0))) + return (x * y) / (x * y); + else + return x; + } + } } - } - } } strong_alias (__ieee754_remainder, __remainder_finite) |