]>
Commit | Line | Data |
---|---|---|
f7eac6eb | 1 | /* |
e4d82761 | 2 | * IBM Accurate Mathematical Library |
aeb25823 | 3 | * written by International Business Machines Corp. |
2b778ceb | 4 | * Copyright (C) 2001-2021 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 | |
5a82c748 | 17 | * along with this program; if not, see <https://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" | |
0e9be4db | 36 | #include <math.h> |
1ed0291c | 37 | #include <math_private.h> |
70e2ba33 | 38 | #include <fenv_private.h> |
220622dd | 39 | #include <libm-alias-finite.h> |
f7eac6eb | 40 | |
e4d82761 UD |
41 | /**************************************************************************/ |
42 | /* An ultimate remainder routine. Given two IEEE double machine numbers x */ | |
43 | /* ,y it computes the correctly rounded (to nearest) value of remainder */ | |
44 | /**************************************************************************/ | |
c5d5d574 OB |
45 | double |
46 | __ieee754_remainder (double x, double y) | |
f7eac6eb | 47 | { |
c5d5d574 OB |
48 | double z, d, xx; |
49 | int4 kx, ky, n, nn, n1, m1, l; | |
50 | mynumber u, t, w = { { 0, 0 } }, v = { { 0, 0 } }, ww = { { 0, 0 } }, r; | |
51 | u.x = x; | |
52 | t.x = y; | |
53 | kx = u.i[HIGH_HALF] & 0x7fffffff; /* no sign for x*/ | |
54 | t.i[HIGH_HALF] &= 0x7fffffff; /*no sign for y */ | |
55 | ky = t.i[HIGH_HALF]; | |
29132b91 | 56 | /*------ |x| < 2^1023 and 2^-970 < |y| < 2^1024 ------------------*/ |
c5d5d574 OB |
57 | if (kx < 0x7fe00000 && ky < 0x7ff00000 && ky >= 0x03500000) |
58 | { | |
59 | SET_RESTORE_ROUND_NOEX (FE_TONEAREST); | |
60 | if (kx + 0x00100000 < ky) | |
61 | return x; | |
62 | if ((kx - 0x01500000) < ky) | |
63 | { | |
64 | z = x / t.x; | |
65 | v.i[HIGH_HALF] = t.i[HIGH_HALF]; | |
66 | d = (z + big.x) - big.x; | |
67 | xx = (x - d * v.x) - d * (t.x - v.x); | |
68 | if (d - z != 0.5 && d - z != -0.5) | |
69 | return (xx != 0) ? xx : ((x > 0) ? ZERO.x : nZERO.x); | |
70 | else | |
71 | { | |
0e9be4db | 72 | if (fabs (xx) > 0.5 * t.x) |
c5d5d574 OB |
73 | return (z > d) ? xx - t.x : xx + t.x; |
74 | else | |
75 | return xx; | |
76 | } | |
77 | } /* (kx<(ky+0x01500000)) */ | |
e4d82761 | 78 | else |
c5d5d574 OB |
79 | { |
80 | r.x = 1.0 / t.x; | |
81 | n = t.i[HIGH_HALF]; | |
82 | nn = (n & 0x7ff00000) + 0x01400000; | |
83 | w.i[HIGH_HALF] = n; | |
84 | ww.x = t.x - w.x; | |
85 | l = (kx - nn) & 0xfff00000; | |
86 | n1 = ww.i[HIGH_HALF]; | |
87 | m1 = r.i[HIGH_HALF]; | |
88 | while (l > 0) | |
89 | { | |
90 | r.i[HIGH_HALF] = m1 - l; | |
91 | z = u.x * r.x; | |
92 | w.i[HIGH_HALF] = n + l; | |
93 | ww.i[HIGH_HALF] = (n1) ? n1 + l : n1; | |
94 | d = (z + big.x) - big.x; | |
95 | u.x = (u.x - d * w.x) - d * ww.x; | |
96 | l = (u.i[HIGH_HALF] & 0x7ff00000) - nn; | |
97 | } | |
98 | r.i[HIGH_HALF] = m1; | |
99 | w.i[HIGH_HALF] = n; | |
100 | ww.i[HIGH_HALF] = n1; | |
101 | z = u.x * r.x; | |
102 | d = (z + big.x) - big.x; | |
103 | u.x = (u.x - d * w.x) - d * ww.x; | |
0e9be4db | 104 | if (fabs (u.x) < 0.5 * t.x) |
c5d5d574 OB |
105 | return (u.x != 0) ? u.x : ((x > 0) ? ZERO.x : nZERO.x); |
106 | else | |
0e9be4db | 107 | if (fabs (u.x) > 0.5 * t.x) |
c5d5d574 OB |
108 | return (d > z) ? u.x + t.x : u.x - t.x; |
109 | else | |
110 | { | |
111 | z = u.x / t.x; d = (z + big.x) - big.x; | |
112 | return ((u.x - d * w.x) - d * ww.x); | |
113 | } | |
114 | } | |
115 | } /* (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000) */ | |
116 | else | |
117 | { | |
118 | if (kx < 0x7fe00000 && ky < 0x7ff00000 && (ky > 0 || t.i[LOW_HALF] != 0)) | |
119 | { | |
0e9be4db | 120 | y = fabs (y) * t128.x; |
c5d5d574 OB |
121 | z = __ieee754_remainder (x, y) * t128.x; |
122 | z = __ieee754_remainder (z, y) * tm128.x; | |
123 | return z; | |
124 | } | |
a1cbf437 | 125 | else |
c5d5d574 OB |
126 | { |
127 | if ((kx & 0x7ff00000) == 0x7fe00000 && ky < 0x7ff00000 && | |
128 | (ky > 0 || t.i[LOW_HALF] != 0)) | |
129 | { | |
0e9be4db | 130 | y = fabs (y); |
c5d5d574 | 131 | z = 2.0 * __ieee754_remainder (0.5 * x, y); |
0e9be4db WD |
132 | d = fabs (z); |
133 | if (d <= fabs (d - y)) | |
c5d5d574 | 134 | return z; |
444ec6b8 JM |
135 | else if (d == y) |
136 | return 0.0 * x; | |
c5d5d574 OB |
137 | else |
138 | return (z > 0) ? z - y : z + y; | |
139 | } | |
140 | else /* if x is too big */ | |
141 | { | |
142 | if (ky == 0 && t.i[LOW_HALF] == 0) /* y = 0 */ | |
143 | return (x * y) / (x * y); | |
144 | else if (kx >= 0x7ff00000 /* x not finite */ | |
145 | || (ky > 0x7ff00000 /* y is NaN */ | |
146 | || (ky == 0x7ff00000 && t.i[LOW_HALF] != 0))) | |
147 | return (x * y) / (x * y); | |
148 | else | |
149 | return x; | |
150 | } | |
151 | } | |
e4d82761 | 152 | } |
f7eac6eb | 153 | } |
220622dd | 154 | libm_alias_finite (__ieee754_remainder, __remainder) |