1 /* Copyright JS Foundation and other contributors, http://js.foundation
2 *
3 * Licensed under the Apache License, Version 2.0 (the "License");
4 * you may not use this file except in compliance with the License.
5 * You may obtain a copy of the License at
6 *
7 * http://www.apache.org/licenses/LICENSE-2.0
8 *
9 * Unless required by applicable law or agreed to in writing, software
10 * distributed under the License is distributed on an "AS IS" BASIS
11 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 * See the License for the specific language governing permissions and
13 * limitations under the License.
14 *
15 * This file is based on work under the following copyright and permission
16 * notice:
17 *
18 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
19 *
20 * Developed at SunSoft, a Sun Microsystems, Inc. business.
21 * Permission to use, copy, modify, and distribute this
22 * software is freely granted, provided that this notice
23 * is preserved.
24 *
25 * @(#)e_fmod.c 1.3 95/01/18
26 */
27
28 #include "jerry-libm-internal.h"
29
30 /* fmod(x,y)
31 * Return x mod y in exact arithmetic
32 *
33 * Method: shift and subtract
34 */
35
36 static const double Zero[] = { 0.0, -0.0, };
37
38 double
fmod(double x,double y)39 fmod (double x, double y)
40 {
41 int n, hx, hy, hz, ix, iy, sx, i;
42 unsigned lx, ly, lz;
43
44 hx = __HI (x); /* high word of x */
45 lx = __LO (x); /* low word of x */
46 hy = __HI (y); /* high word of y */
47 ly = __LO (y); /* low word of y */
48 sx = hx & 0x80000000; /* sign of x */
49 hx ^= sx; /* |x| */
50 hy &= 0x7fffffff; /* |y| */
51
52 /* purge off exception values */
53 if ((hy | ly) == 0 || (hx >= 0x7ff00000) || /* y = 0, or x not finite */
54 ((hy | ((ly | -ly) >> 31)) > 0x7ff00000)) /* or y is NaN */
55 {
56 return NAN;
57 }
58 if (hx <= hy)
59 {
60 if ((hx < hy) || (lx < ly)) /* |x| < |y| return x */
61 {
62 return x;
63 }
64 if (lx == ly) /* |x| = |y| return x * 0 */
65 {
66 return Zero[(unsigned) sx >> 31];
67 }
68 }
69
70 /* determine ix = ilogb(x) */
71 if (hx < 0x00100000) /* subnormal x */
72 {
73 if (hx == 0)
74 {
75 for (ix = -1043, i = lx; i > 0; i <<= 1)
76 {
77 ix -= 1;
78 }
79 }
80 else
81 {
82 for (ix = -1022, i = (hx << 11); i > 0; i <<= 1)
83 {
84 ix -= 1;
85 }
86 }
87 }
88 else
89 {
90 ix = (hx >> 20) - 1023;
91 }
92
93 /* determine iy = ilogb(y) */
94 if (hy < 0x00100000) /* subnormal y */
95 {
96 if (hy == 0)
97 {
98 for (iy = -1043, i = ly; i > 0; i <<= 1)
99 {
100 iy -= 1;
101 }
102 }
103 else
104 {
105 for (iy = -1022, i = (hy << 11); i > 0; i <<= 1)
106 {
107 iy -= 1;
108 }
109 }
110 }
111 else
112 {
113 iy = (hy >> 20) - 1023;
114 }
115
116 /* set up {hx,lx}, {hy,ly} and align y to x */
117 if (ix >= -1022)
118 {
119 hx = 0x00100000 | (0x000fffff & hx);
120 }
121 else /* subnormal x, shift x to normal */
122 {
123 n = -1022 - ix;
124 if (n <= 31)
125 {
126 hx = (((unsigned int) hx) << n) | (lx >> (32 - n));
127 lx <<= n;
128 }
129 else
130 {
131 hx = lx << (n - 32);
132 lx = 0;
133 }
134 }
135 if (iy >= -1022)
136 {
137 hy = 0x00100000 | (0x000fffff & hy);
138 }
139 else /* subnormal y, shift y to normal */
140 {
141 n = -1022 - iy;
142 if (n <= 31)
143 {
144 hy = (((unsigned int) hy) << n) | (ly >> (32 - n));
145 ly <<= n;
146 }
147 else
148 {
149 hy = ly << (n - 32);
150 ly = 0;
151 }
152 }
153
154 /* fix point fmod */
155 n = ix - iy;
156 while (n--)
157 {
158 hz = hx - hy;
159 lz = lx - ly;
160 if (lx < ly)
161 {
162 hz -= 1;
163 }
164 if (hz < 0)
165 {
166 hx = hx + hx + (lx >> 31);
167 lx = lx + lx;
168 }
169 else
170 {
171 if ((hz | lz) == 0) /* return sign(x) * 0 */
172 {
173 return Zero[(unsigned) sx >> 31];
174 }
175 hx = hz + hz + (lz >> 31);
176 lx = lz + lz;
177 }
178 }
179 hz = hx - hy;
180 lz = lx - ly;
181 if (lx < ly)
182 {
183 hz -= 1;
184 }
185 if (hz >= 0)
186 {
187 hx = hz;
188 lx = lz;
189 }
190
191 /* convert back to floating value and restore the sign */
192 if ((hx | lx) == 0) /* return sign(x) * 0 */
193 {
194 return Zero[(unsigned) sx >> 31];
195 }
196 while (hx < 0x00100000) /* normalize x */
197 {
198 hx = hx + hx + (lx >> 31);
199 lx = lx + lx;
200 iy -= 1;
201 }
202
203 double_accessor ret;
204 if (iy >= -1022) /* normalize output */
205 {
206 hx = ((hx - 0x00100000) | ((iy + 1023) << 20));
207 ret.as_int.hi = hx | sx;
208 ret.as_int.lo = lx;
209 }
210 else /* subnormal output */
211 {
212 n = -1022 - iy;
213 if (n <= 20)
214 {
215 lx = (lx >> n) | ((unsigned) hx << (32 - n));
216 hx >>= n;
217 }
218 else if (n <= 31)
219 {
220 lx = (hx << (32 - n)) | (lx >> n);
221 hx = sx;
222 }
223 else
224 {
225 lx = hx >> (n - 32);
226 hx = sx;
227 }
228 ret.as_int.hi = hx | sx;
229 ret.as_int.lo = lx;
230 }
231 return ret.dbl; /* exact output */
232 } /* fmod */
233