This source file includes following definitions.
- FMOD
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 #if ! defined USE_LONG_DOUBLE
18 # include <config.h>
19 #endif
20
21
22 #include <math.h>
23
24 #include <float.h>
25 #include <stdlib.h>
26
27 #ifdef USE_LONG_DOUBLE
28 # define FMOD fmodl
29 # define DOUBLE long double
30 # define MANT_DIG LDBL_MANT_DIG
31 # define L_(literal) literal##L
32 # define FREXP frexpl
33 # define LDEXP ldexpl
34 # define FABS fabsl
35 # define TRUNC truncl
36 # define ISNAN isnanl
37 #else
38 # define FMOD fmod
39 # define DOUBLE double
40 # define MANT_DIG DBL_MANT_DIG
41 # define L_(literal) literal
42 # define FREXP frexp
43 # define LDEXP ldexp
44 # define FABS fabs
45 # define TRUNC trunc
46 # define ISNAN isnand
47 #endif
48
49
50
51 #if defined _MSC_VER && !defined __clang__
52 # pragma fenv_access (off)
53 #endif
54
55 #undef NAN
56 #if defined _MSC_VER
57 static DOUBLE zero;
58 # define NAN (zero / zero)
59 #else
60 # define NAN (L_(0.0) / L_(0.0))
61 #endif
62
63
64
65
66
67
68
69
70
71
72
73
74
75 #define LIMB_BITS (MANT_DIG / 2)
76
77
78 static const DOUBLE TWO_LIMB_BITS =
79
80
81
82 (DOUBLE) (1U << (LIMB_BITS / 3))
83 * (DOUBLE) (1U << ((LIMB_BITS + 1) / 3))
84 * (DOUBLE) (1U << ((LIMB_BITS + 2) / 3));
85
86
87 static const DOUBLE TWO_LIMB_BITS_INVERSE =
88
89
90
91 L_(1.0) / ((DOUBLE) (1U << (LIMB_BITS / 3))
92 * (DOUBLE) (1U << ((LIMB_BITS + 1) / 3))
93 * (DOUBLE) (1U << ((LIMB_BITS + 2) / 3)));
94
95 DOUBLE
96 FMOD (DOUBLE x, DOUBLE y)
97 {
98 if (isfinite (x) && isfinite (y) && y != L_(0.0))
99 {
100 if (x == L_(0.0))
101
102 return x;
103
104 {
105 int negate = ((!signbit (x)) ^ (!signbit (y)));
106
107
108 x = FABS (x);
109 y = FABS (y);
110
111
112 if (x < y)
113 return (negate ? - x : x);
114
115 {
116 int yexp;
117 DOUBLE ym;
118 DOUBLE y1;
119 DOUBLE y0;
120 int k;
121 DOUBLE x2;
122 DOUBLE x1;
123 DOUBLE x0;
124
125
126
127
128
129
130 ym = FREXP (y, &yexp);
131 ym = ym * TWO_LIMB_BITS;
132 y1 = TRUNC (ym);
133 y0 = (ym - y1) * TWO_LIMB_BITS;
134
135
136
137
138
139 {
140 int xexp;
141 DOUBLE xm = FREXP (x, &xexp);
142
143 xexp -= yexp;
144
145 k = (xexp + LIMB_BITS - 1) / LIMB_BITS;
146
147
148 xm = LDEXP (xm, xexp - (k - 1) * LIMB_BITS);
149
150
151
152 x2 = TRUNC (xm);
153 x1 = (xm - x2) * TWO_LIMB_BITS;
154
155 }
156
157
158 if (x2 > y1 || (x2 == y1 && x1 >= y0))
159 {
160
161 x2 -= y1;
162 x1 -= y0;
163 if (x1 < L_(0.0))
164 {
165 if (!(x2 >= L_(1.0)))
166 abort ();
167 x2 -= L_(1.0);
168 x1 += TWO_LIMB_BITS;
169 }
170 }
171
172
173 {
174 DOUBLE x1int = TRUNC (x1);
175 x0 = TRUNC ((x1 - x1int) * TWO_LIMB_BITS);
176 x1 = x1int;
177 }
178
179 for (; k > 0; k--)
180 {
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224 DOUBLE q =
225 (x2 >= y1
226 ? TWO_LIMB_BITS - L_(1.0)
227 : TRUNC ((x2 * TWO_LIMB_BITS + x1) / y1));
228 if (q > L_(0.0))
229 {
230
231
232
233 DOUBLE q_y1 = q * y1;
234 DOUBLE q_y1_1 = TRUNC (q_y1 * TWO_LIMB_BITS_INVERSE);
235 DOUBLE q_y1_0 = q_y1 - q_y1_1 * TWO_LIMB_BITS;
236 DOUBLE q_y0 = q * y0;
237 DOUBLE q_y0_1 = TRUNC (q_y0 * TWO_LIMB_BITS_INVERSE);
238 DOUBLE q_y0_0 = q_y0 - q_y0_1 * TWO_LIMB_BITS;
239 x2 -= q_y1_1;
240 x1 -= q_y1_0;
241 x1 -= q_y0_1;
242 x0 -= q_y0_0;
243
244 if (x0 < L_(0.0))
245 {
246 x0 += TWO_LIMB_BITS;
247 x1 -= L_(1.0);
248 }
249 if (x1 < L_(0.0))
250 {
251 x1 += TWO_LIMB_BITS;
252 x2 -= L_(1.0);
253 if (x1 < L_(0.0))
254 {
255 x1 += TWO_LIMB_BITS;
256 x2 -= L_(1.0);
257 }
258 }
259 if (x2 < L_(0.0))
260 {
261
262 x1 += y1;
263 x0 += y0;
264
265 if (x0 >= TWO_LIMB_BITS)
266 {
267 x0 -= TWO_LIMB_BITS;
268 x1 += L_(1.0);
269 }
270 if (x1 >= TWO_LIMB_BITS)
271 {
272 x1 -= TWO_LIMB_BITS;
273 x2 += L_(1.0);
274 }
275 if (x2 < L_(0.0))
276 {
277
278 x1 += y1;
279 x0 += y0;
280
281 if (x0 >= TWO_LIMB_BITS)
282 {
283 x0 -= TWO_LIMB_BITS;
284 x1 += L_(1.0);
285 }
286 if (x1 >= TWO_LIMB_BITS)
287 {
288 x1 -= TWO_LIMB_BITS;
289 x2 += L_(1.0);
290 }
291 if (x2 < L_(0.0))
292
293
294 abort ();
295 }
296 }
297 }
298 if (x2 > L_(0.0)
299 || x1 > y1
300 || (x1 == y1 && x0 >= y0))
301 {
302
303 x1 -= y1;
304 x0 -= y0;
305
306 if (x0 < L_(0.0))
307 {
308 x0 += TWO_LIMB_BITS;
309 x1 -= L_(1.0);
310 }
311 if (x1 < L_(0.0))
312 {
313 x1 += TWO_LIMB_BITS;
314 x2 -= L_(1.0);
315 }
316 if (x2 < L_(0.0))
317 abort ();
318 if (x2 > L_(0.0)
319 || x1 > y1
320 || (x1 == y1 && x0 >= y0))
321
322
323 abort ();
324 }
325
326
327 x2 = x1;
328 #if (MANT_DIG + 1) / 2 > LIMB_BITS
329 x1 = TRUNC (x0);
330 x0 = (x0 - x1) * TWO_LIMB_BITS;
331 #else
332 x1 = x0;
333 x0 = L_(0.0);
334 #endif
335
336 }
337
338
339
340
341 {
342 DOUBLE r =
343 LDEXP ((x2 * TWO_LIMB_BITS + x1) * TWO_LIMB_BITS + x0,
344 yexp - 3 * LIMB_BITS);
345 return (negate ? - r : r);
346 }
347 }
348 }
349 }
350 else
351 {
352 if (ISNAN (x) || ISNAN (y))
353 return x + y;
354 else if (isinf (y))
355 return x;
356 else
357
358 return NAN;
359 }
360 }