1 | /* |
2 | * IBM Accurate Mathematical Library |
3 | * written by International Business Machines Corp. |
4 | * Copyright (C) 2001-2018 Free Software Foundation, Inc. |
5 | * |
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 |
8 | * the Free Software Foundation; either version 2.1 of the License, or |
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 |
14 | * GNU Lesser General Public License for more details. |
15 | * |
16 | * You should have received a copy of the GNU Lesser General Public License |
17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
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.h> |
37 | #include <math_private.h> |
38 | |
39 | /**************************************************************************/ |
40 | /* An ultimate remainder routine. Given two IEEE double machine numbers x */ |
41 | /* ,y it computes the correctly rounded (to nearest) value of remainder */ |
42 | /**************************************************************************/ |
43 | double |
44 | __ieee754_remainder (double x, double y) |
45 | { |
46 | double z, d, xx; |
47 | int4 kx, ky, n, nn, n1, m1, l; |
48 | mynumber u, t, w = { { 0, 0 } }, v = { { 0, 0 } }, ww = { { 0, 0 } }, r; |
49 | u.x = x; |
50 | t.x = y; |
51 | kx = u.i[HIGH_HALF] & 0x7fffffff; /* no sign for x*/ |
52 | t.i[HIGH_HALF] &= 0x7fffffff; /*no sign for y */ |
53 | ky = t.i[HIGH_HALF]; |
54 | /*------ |x| < 2^1023 and 2^-970 < |y| < 2^1024 ------------------*/ |
55 | if (kx < 0x7fe00000 && ky < 0x7ff00000 && ky >= 0x03500000) |
56 | { |
57 | SET_RESTORE_ROUND_NOEX (FE_TONEAREST); |
58 | if (kx + 0x00100000 < ky) |
59 | return x; |
60 | if ((kx - 0x01500000) < ky) |
61 | { |
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) |
67 | return (xx != 0) ? xx : ((x > 0) ? ZERO.x : nZERO.x); |
68 | else |
69 | { |
70 | if (fabs (xx) > 0.5 * t.x) |
71 | return (z > d) ? xx - t.x : xx + t.x; |
72 | else |
73 | return xx; |
74 | } |
75 | } /* (kx<(ky+0x01500000)) */ |
76 | else |
77 | { |
78 | r.x = 1.0 / t.x; |
79 | n = t.i[HIGH_HALF]; |
80 | nn = (n & 0x7ff00000) + 0x01400000; |
81 | w.i[HIGH_HALF] = n; |
82 | ww.x = t.x - w.x; |
83 | l = (kx - nn) & 0xfff00000; |
84 | n1 = ww.i[HIGH_HALF]; |
85 | m1 = r.i[HIGH_HALF]; |
86 | while (l > 0) |
87 | { |
88 | r.i[HIGH_HALF] = m1 - l; |
89 | z = u.x * r.x; |
90 | w.i[HIGH_HALF] = n + l; |
91 | ww.i[HIGH_HALF] = (n1) ? n1 + l : n1; |
92 | d = (z + big.x) - big.x; |
93 | u.x = (u.x - d * w.x) - d * ww.x; |
94 | l = (u.i[HIGH_HALF] & 0x7ff00000) - nn; |
95 | } |
96 | r.i[HIGH_HALF] = m1; |
97 | w.i[HIGH_HALF] = n; |
98 | ww.i[HIGH_HALF] = n1; |
99 | z = u.x * r.x; |
100 | d = (z + big.x) - big.x; |
101 | u.x = (u.x - d * w.x) - d * ww.x; |
102 | if (fabs (u.x) < 0.5 * t.x) |
103 | return (u.x != 0) ? u.x : ((x > 0) ? ZERO.x : nZERO.x); |
104 | else |
105 | if (fabs (u.x) > 0.5 * t.x) |
106 | return (d > z) ? u.x + t.x : u.x - t.x; |
107 | else |
108 | { |
109 | z = u.x / t.x; d = (z + big.x) - big.x; |
110 | return ((u.x - d * w.x) - d * ww.x); |
111 | } |
112 | } |
113 | } /* (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000) */ |
114 | else |
115 | { |
116 | if (kx < 0x7fe00000 && ky < 0x7ff00000 && (ky > 0 || t.i[LOW_HALF] != 0)) |
117 | { |
118 | y = fabs (y) * t128.x; |
119 | z = __ieee754_remainder (x, y) * t128.x; |
120 | z = __ieee754_remainder (z, y) * tm128.x; |
121 | return z; |
122 | } |
123 | else |
124 | { |
125 | if ((kx & 0x7ff00000) == 0x7fe00000 && ky < 0x7ff00000 && |
126 | (ky > 0 || t.i[LOW_HALF] != 0)) |
127 | { |
128 | y = fabs (y); |
129 | z = 2.0 * __ieee754_remainder (0.5 * x, y); |
130 | d = fabs (z); |
131 | if (d <= fabs (d - y)) |
132 | return z; |
133 | else if (d == y) |
134 | return 0.0 * x; |
135 | else |
136 | return (z > 0) ? z - y : z + y; |
137 | } |
138 | else /* if x is too big */ |
139 | { |
140 | if (ky == 0 && t.i[LOW_HALF] == 0) /* y = 0 */ |
141 | return (x * y) / (x * y); |
142 | else if (kx >= 0x7ff00000 /* x not finite */ |
143 | || (ky > 0x7ff00000 /* y is NaN */ |
144 | || (ky == 0x7ff00000 && t.i[LOW_HALF] != 0))) |
145 | return (x * y) / (x * y); |
146 | else |
147 | return x; |
148 | } |
149 | } |
150 | } |
151 | } |
152 | strong_alias (__ieee754_remainder, __remainder_finite) |
153 | |