2 * IBM Accurate Mathematical Library
3 * Copyright (c) International Business Machines Corp., 2001
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU Lesser General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
15 * You should have received a copy of the GNU Lesser General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
19 /**************************************************************************/
20 /* MODULE_NAME urem.c */
22 /* FUNCTION: uremainder */
24 /* An ultimate remainder routine. Given two IEEE double machine numbers x */
25 /* ,y it computes the correctly rounded (to nearest) value of remainder */
26 /* of dividing x by y. */
27 /* Assumption: Machine arithmetic operations are performed in */
28 /* round to nearest mode of IEEE 754 standard. */
30 /* ************************************************************************/
36 #include "math_private.h"
38 /**************************************************************************/
39 /* An ultimate remainder routine. Given two IEEE double machine numbers x */
40 /* ,y it computes the correctly rounded (to nearest) value of remainder */
41 /**************************************************************************/
42 double __ieee754_remainder(double x, double y)
48 int4 kx,ky,n,nn,n1,m1,l;
52 mynumber u,t,w={{0,0}},v={{0,0}},ww={{0,0}},r;
55 kx=u.i[HIGH_HALF]&0x7fffffff; /* no sign for x*/
56 t.i[HIGH_HALF]&=0x7fffffff; /*no sign for y */
58 /*------ |x| < 2^1023 and 2^-970 < |y| < 2^1024 ------------------*/
59 if (kx<0x7fe00000 && ky<0x7ff00000 && ky>=0x03500000) {
60 if (kx+0x00100000<ky) return x;
61 if ((kx-0x01500000)<ky) {
63 v.i[HIGH_HALF]=t.i[HIGH_HALF];
65 xx=(x-d*v.x)-d*(t.x-v.x);
66 if (d-z!=0.5&&d-z!=-0.5) return (xx!=0)?xx:((x>0)?ZERO.x:nZERO.x);
68 if (ABS(xx)>0.5*t.x) return (z>d)?xx-t.x:xx+t.x;
71 } /* (kx<(ky+0x01500000)) */
75 nn=(n&0x7ff00000)+0x01400000;
85 ww.i[HIGH_HALF]=(n1)?n1+l:n1;
87 u.x=(u.x-d*w.x)-d*ww.x;
88 l=(u.i[HIGH_HALF]&0x7ff00000)-nn;
95 u.x=(u.x-d*w.x)-d*ww.x;
96 if (ABS(u.x)<0.5*t.x) return (u.x!=0)?u.x:((x>0)?ZERO.x:nZERO.x);
98 if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x;
100 {z=u.x/t.x; d=(z+big.x)-big.x; return ((u.x-d*w.x)-d*ww.x);}
103 } /* (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000) */
105 if (kx<0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
107 z=__ieee754_remainder(x,y)*t128.x;
108 z=__ieee754_remainder(z,y)*tm128.x;
112 if ((kx&0x7ff00000)==0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
114 z=2.0*__ieee754_remainder(0.5*x,y);
116 if (d <= ABS(d-y)) return z;
117 else return (z>0)?z-y:z+y;
119 else { /* if x is too big */
120 if (kx == 0x7ff00000 && u.i[LOW_HALF] == 0 && y == 1.0)
122 if (kx>=0x7ff00000||(ky==0&&t.i[LOW_HALF]==0)||ky>0x7ff00000||
123 (ky==0x7ff00000&&t.i[LOW_HALF]!=0))
124 return (u.i[HIGH_HALF]&0x80000000)?nNAN.x:NAN.x;