]>
Commit | Line | Data |
---|---|---|
f7eac6eb | 1 | /* |
e4d82761 | 2 | * IBM Accurate Mathematical Library |
aeb25823 | 3 | * written by International Business Machines Corp. |
568035b7 | 4 | * Copyright (C) 2001-2013 Free Software Foundation, Inc. |
0ac5ae23 | 5 | * |
e4d82761 UD |
6 | * This program is free software; you can redistribute it and/or modify |
7 | * it under the terms of the GNU Lesser General Public License as published by | |
cc7375ce | 8 | * the Free Software Foundation; either version 2.1 of the License, or |
e4d82761 UD |
9 | * (at your option) any later version. |
10 | * | |
11 | * This program is distributed in the hope that it will be useful, | |
12 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
c6c6dd48 | 14 | * GNU Lesser General Public License for more details. |
e4d82761 UD |
15 | * |
16 | * You should have received a copy of the GNU Lesser General Public License | |
59ba27a6 | 17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
f7eac6eb | 18 | */ |
e4d82761 UD |
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 | /* ************************************************************************/ | |
f7eac6eb | 31 | |
e4d82761 UD |
32 | #include "endian.h" |
33 | #include "mydefs.h" | |
34 | #include "urem.h" | |
35 | #include "MathLib.h" | |
1ed0291c | 36 | #include <math_private.h> |
f7eac6eb | 37 | |
e4d82761 UD |
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) | |
f7eac6eb | 43 | { |
50944bca | 44 | double z,d,xx; |
50944bca | 45 | int4 kx,ky,n,nn,n1,m1,l; |
50944bca | 46 | mynumber u,t,w={{0,0}},v={{0,0}},ww={{0,0}},r; |
e4d82761 UD |
47 | u.x=x; |
48 | t.x=y; | |
49 | kx=u.i[HIGH_HALF]&0x7fffffff; /* no sign for x*/ | |
50 | t.i[HIGH_HALF]&=0x7fffffff; /*no sign for y */ | |
51 | ky=t.i[HIGH_HALF]; | |
29132b91 UD |
52 | /*------ |x| < 2^1023 and 2^-970 < |y| < 2^1024 ------------------*/ |
53 | if (kx<0x7fe00000 && ky<0x7ff00000 && ky>=0x03500000) { | |
e4d82761 UD |
54 | if (kx+0x00100000<ky) return x; |
55 | if ((kx-0x01500000)<ky) { | |
56 | z=x/t.x; | |
57 | v.i[HIGH_HALF]=t.i[HIGH_HALF]; | |
58 | d=(z+big.x)-big.x; | |
59 | xx=(x-d*v.x)-d*(t.x-v.x); | |
60 | if (d-z!=0.5&&d-z!=-0.5) return (xx!=0)?xx:((x>0)?ZERO.x:nZERO.x); | |
61 | else { | |
62 | if (ABS(xx)>0.5*t.x) return (z>d)?xx-t.x:xx+t.x; | |
63 | else return xx; | |
64 | } | |
65 | } /* (kx<(ky+0x01500000)) */ | |
66 | else { | |
67 | r.x=1.0/t.x; | |
68 | n=t.i[HIGH_HALF]; | |
69 | nn=(n&0x7ff00000)+0x01400000; | |
70 | w.i[HIGH_HALF]=n; | |
71 | ww.x=t.x-w.x; | |
72 | l=(kx-nn)&0xfff00000; | |
73 | n1=ww.i[HIGH_HALF]; | |
74 | m1=r.i[HIGH_HALF]; | |
75 | while (l>0) { | |
76 | r.i[HIGH_HALF]=m1-l; | |
77 | z=u.x*r.x; | |
78 | w.i[HIGH_HALF]=n+l; | |
79 | ww.i[HIGH_HALF]=(n1)?n1+l:n1; | |
80 | d=(z+big.x)-big.x; | |
81 | u.x=(u.x-d*w.x)-d*ww.x; | |
82 | l=(u.i[HIGH_HALF]&0x7ff00000)-nn; | |
83 | } | |
84 | r.i[HIGH_HALF]=m1; | |
85 | w.i[HIGH_HALF]=n; | |
86 | ww.i[HIGH_HALF]=n1; | |
87 | z=u.x*r.x; | |
88 | d=(z+big.x)-big.x; | |
89 | u.x=(u.x-d*w.x)-d*ww.x; | |
90 | if (ABS(u.x)<0.5*t.x) return (u.x!=0)?u.x:((x>0)?ZERO.x:nZERO.x); | |
91 | else | |
0ac5ae23 UD |
92 | if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x; |
93 | else | |
94 | {z=u.x/t.x; d=(z+big.x)-big.x; return ((u.x-d*w.x)-d*ww.x);} | |
e4d82761 | 95 | } |
f7eac6eb | 96 | |
29132b91 | 97 | } /* (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000) */ |
e4d82761 | 98 | else { |
29132b91 | 99 | if (kx<0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) { |
e4d82761 | 100 | y=ABS(y)*t128.x; |
ca58f1db UD |
101 | z=__ieee754_remainder(x,y)*t128.x; |
102 | z=__ieee754_remainder(z,y)*tm128.x; | |
e4d82761 UD |
103 | return z; |
104 | } | |
29132b91 UD |
105 | else { |
106 | if ((kx&0x7ff00000)==0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) { | |
107 | y=ABS(y); | |
108 | z=2.0*__ieee754_remainder(0.5*x,y); | |
109 | d = ABS(z); | |
110 | if (d <= ABS(d-y)) return z; | |
111 | else return (z>0)?z-y:z+y; | |
112 | } | |
e4d82761 | 113 | else { /* if x is too big */ |
a1cbf437 TS |
114 | if (ky==0 && t.i[LOW_HALF] == 0) /* y = 0 */ |
115 | return (x * y) / (x * y); | |
116 | else if (kx >= 0x7ff00000 /* x not finite */ | |
117 | || (ky>0x7ff00000 /* y is NaN */ | |
118 | || (ky == 0x7ff00000 && t.i[LOW_HALF] != 0))) | |
119 | return (x * y) / (x * y); | |
120 | else | |
121 | return x; | |
e4d82761 | 122 | } |
29132b91 | 123 | } |
e4d82761 | 124 | } |
f7eac6eb | 125 | } |
0ac5ae23 | 126 | strong_alias (__ieee754_remainder, __remainder_finite) |