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