(__ieee754_remainder): Fix overflow problem.
[kopensolaris-gnu/glibc.git] / sysdeps / ieee754 / dbl-64 / e_remainder.c
1 /*
2  * IBM Accurate Mathematical Library
3  * Copyright (c) International Business Machines Corp., 2001
4  *
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.
9  *
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.
14  *
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.
18  */
19 /**************************************************************************/
20 /*  MODULE_NAME urem.c                                                    */
21 /*                                                                        */
22 /*  FUNCTION: uremainder                                                  */
23 /*                                                                        */
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.                            */
29 /*                                                                        */
30 /* ************************************************************************/
31
32 #include "endian.h"
33 #include "mydefs.h"
34 #include "urem.h"
35 #include "MathLib.h"
36 #include "math_private.h"
37
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)
43 {
44   double z,d,xx;
45 #if 0
46   double yy;
47 #endif
48   int4 kx,ky,n,nn,n1,m1,l;
49 #if 0
50   int4 m;
51 #endif
52   mynumber u,t,w={{0,0}},v={{0,0}},ww={{0,0}},r;
53   u.x=x;
54   t.x=y;
55   kx=u.i[HIGH_HALF]&0x7fffffff; /* no sign  for x*/
56   t.i[HIGH_HALF]&=0x7fffffff;   /*no sign for y */
57   ky=t.i[HIGH_HALF];
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) {
62       z=x/t.x;
63       v.i[HIGH_HALF]=t.i[HIGH_HALF];
64       d=(z+big.x)-big.x;
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);
67       else {
68         if (ABS(xx)>0.5*t.x) return (z>d)?xx-t.x:xx+t.x;
69         else return xx;
70       }
71     }   /*    (kx<(ky+0x01500000))         */
72     else  {
73       r.x=1.0/t.x;
74       n=t.i[HIGH_HALF];
75       nn=(n&0x7ff00000)+0x01400000;
76       w.i[HIGH_HALF]=n;
77       ww.x=t.x-w.x;
78       l=(kx-nn)&0xfff00000;
79       n1=ww.i[HIGH_HALF];
80       m1=r.i[HIGH_HALF];
81       while (l>0) {
82         r.i[HIGH_HALF]=m1-l;
83         z=u.x*r.x;
84         w.i[HIGH_HALF]=n+l;
85         ww.i[HIGH_HALF]=(n1)?n1+l:n1;
86         d=(z+big.x)-big.x;
87         u.x=(u.x-d*w.x)-d*ww.x;
88         l=(u.i[HIGH_HALF]&0x7ff00000)-nn;
89       }
90       r.i[HIGH_HALF]=m1;
91       w.i[HIGH_HALF]=n;
92       ww.i[HIGH_HALF]=n1;
93       z=u.x*r.x;
94       d=(z+big.x)-big.x;
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);
97       else
98         if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x;
99         else
100         {z=u.x/t.x; d=(z+big.x)-big.x; return ((u.x-d*w.x)-d*ww.x);}
101     }
102
103   }   /*   (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000)     */
104   else {
105     if (kx<0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
106       y=ABS(y)*t128.x;
107       z=__ieee754_remainder(x,y)*t128.x;
108       z=__ieee754_remainder(z,y)*tm128.x;
109       return z;
110     }
111   else {
112     if ((kx&0x7ff00000)==0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
113       y=ABS(y);
114       z=2.0*__ieee754_remainder(0.5*x,y);
115       d = ABS(z);
116       if (d <= ABS(d-y)) return z;
117       else return (z>0)?z-y:z+y;
118     }
119     else { /* if x is too big */
120       if (kx == 0x7ff00000 && u.i[LOW_HALF] == 0 && y == 1.0)
121         return x / x;
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;
125       else return x;
126     }
127    }
128   }
129 }