1 /* tfmod -- test file for mpfr_fmod
2
3 Copyright 2007-2020 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 #include "mpfr-test.h"
24
25 #define TEST_FUNCTION mpfr_fmod
26 #define TWO_ARGS
27 #include "tgeneric.c"
28
29 /* compute remainder as in definition:
30 r = x - n * y, where n = trunc(x/y).
31 warning: may change flags. */
32 static int
slow_fmod(mpfr_ptr r,mpfr_srcptr x,mpfr_srcptr y,mpfr_rnd_t rnd)33 slow_fmod (mpfr_ptr r, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
34 {
35 mpfr_t q;
36 int inexact;
37 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x) || MPFR_IS_SINGULAR (y)))
38 {
39 if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y) || MPFR_IS_INF (x)
40 || MPFR_IS_ZERO (y))
41 {
42 MPFR_SET_NAN (r);
43 MPFR_RET_NAN;
44 }
45 else /* either y is Inf and x is 0 or non-special,
46 or x is 0 and y is non-special,
47 in both cases the quotient is zero. */
48 return mpfr_set (r, x, rnd);
49 }
50 /* regular cases */
51 /* if 2^(ex-1) <= |x| < 2^ex, and 2^(ey-1) <= |y| < 2^ey,
52 then |x/y| < 2^(ex-ey+1) */
53 mpfr_init2 (q,
54 MAX (MPFR_PREC_MIN, mpfr_get_exp (x) - mpfr_get_exp (y) + 1));
55 mpfr_div (q, x, y, MPFR_RNDZ);
56 mpfr_trunc (q, q); /* may change inexact flag */
57 mpfr_prec_round (q, mpfr_get_prec (q) + mpfr_get_prec (y), MPFR_RNDZ);
58 inexact = mpfr_mul (q, q, y, MPFR_RNDZ); /* exact */
59 inexact = mpfr_sub (r, x, q, rnd);
60 mpfr_clear (q);
61 return inexact;
62 }
63
64 static void
test_failed(mpfr_ptr erem,mpfr_ptr grem,int eret,int gret,mpfr_ptr x,mpfr_ptr y,mpfr_rnd_t rnd)65 test_failed (mpfr_ptr erem, mpfr_ptr grem, int eret, int gret,
66 mpfr_ptr x, mpfr_ptr y, mpfr_rnd_t rnd)
67 {
68 printf ("error: mpfr_fmod (r, x, y, rnd)\n x = ");
69 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDD);
70 printf ("\n y = ");
71 mpfr_out_str (stdout, 10, 0, y, MPFR_RNDD);
72 printf ("\nrnd = %s", mpfr_print_rnd_mode (rnd));
73 if (eret != gret)
74 printf ("\nexpected %s return value, got %d",
75 (eret < 0 ? "negative" : eret > 0 ? "positive" : "zero"), gret);
76 printf ("\n expected r = ");
77 mpfr_out_str (stdout, 10, 0, erem, MPFR_RNDD);
78 printf ("\n got r = ");
79 mpfr_out_str (stdout, 10, 0, grem, MPFR_RNDD);
80 putchar ('\n');
81
82 exit (1);
83 }
84
85 static void
check(mpfr_ptr r0,mpfr_ptr x,mpfr_ptr y,mpfr_rnd_t rnd)86 check (mpfr_ptr r0, mpfr_ptr x, mpfr_ptr y, mpfr_rnd_t rnd)
87 {
88 int inex0, inex1;
89 mpfr_t r1;
90 mpfr_init2 (r1, mpfr_get_prec (r0));
91
92 inex0 = mpfr_fmod (r0, x, y, rnd);
93 inex1 = slow_fmod (r1, x, y, rnd);
94 if (!mpfr_equal_p (r0, r1) || inex0 != inex1)
95 test_failed (r1, r0, inex1, inex0, x, y, rnd);
96 mpfr_clear (r1);
97 }
98
99 static void
regular(void)100 regular (void)
101 {
102 mpfr_t x, y, r;
103 mpfr_inits (x, y, r, (mpfr_ptr) 0);
104
105 /* remainder = 0 */
106 mpfr_set_str (y, "FEDCBA987654321p-64", 16, MPFR_RNDN);
107 mpfr_pow_ui (x, y, 42, MPFR_RNDN);
108 check (r, x, y, MPFR_RNDN);
109
110 /* x < y */
111 mpfr_set_ui_2exp (x, 64723, -19, MPFR_RNDN);
112 mpfr_mul (x, x, y, MPFR_RNDN);
113 check (r, x, y, MPFR_RNDN);
114
115 /* sign(x) = sign (r) */
116 mpfr_set_ui (x, 123798, MPFR_RNDN);
117 mpfr_set_ui (y, 10, MPFR_RNDN);
118 check (r, x, y, MPFR_RNDN);
119
120 /* huge difference between precisions */
121 mpfr_set_prec (x, 314);
122 mpfr_set_prec (y, 8);
123 mpfr_set_prec (r, 123);
124 mpfr_const_pi (x, MPFR_RNDD); /* x = pi */
125 mpfr_set_ui_2exp (y, 1, 3, MPFR_RNDD); /* y = 1/8 */
126 check (r, x, y, MPFR_RNDD);
127
128 mpfr_clears (x, y, r, (mpfr_ptr) 0);
129 }
130
131 static void
special(void)132 special (void)
133 {
134 int inexact;
135 mpfr_t x, y, r, t;
136
137 mpfr_inits (x, y, r, t, (mpfr_ptr) 0);
138
139 mpfr_set_nan (t);
140
141 /* fmod (NaN, NaN) is NaN */
142 mpfr_set_nan (x);
143 mpfr_set_nan (y);
144 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
145 if (!mpfr_nan_p (r) || inexact != 0)
146 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN);
147
148 /* fmod (NaN, +0) is NaN */
149 mpfr_set_ui (y, 0, MPFR_RNDN);
150 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
151 if (!mpfr_nan_p (r) || inexact != 0)
152 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN);
153
154 /* fmod (+1, 0) is NaN */
155 mpfr_set_ui (x, 1, MPFR_RNDN);
156 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
157 if (!mpfr_nan_p (r) || inexact != 0)
158 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN);
159
160 /* fmod (0, 0) is NaN */
161 mpfr_set_ui (x, 0, MPFR_RNDN);
162 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
163 if (!mpfr_nan_p (r) || inexact != 0)
164 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN);
165
166 /* fmod (+inf, +0) is NaN */
167 mpfr_set_inf (x, +1);
168 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
169 if (!mpfr_nan_p (r) || inexact != 0)
170 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN);
171
172 /* fmod (-inf, +0) is NaN */
173 mpfr_set_inf (x, -1);
174 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
175 if (!mpfr_nan_p (r) || inexact != 0)
176 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN);
177
178 /* fmod (-inf, -0) is NaN */
179 mpfr_neg (x, x, MPFR_RNDN);
180 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
181 if (!mpfr_nan_p (r) || inexact != 0)
182 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN);
183
184 /* fmod (-inf, +1) is NaN */
185 mpfr_set_ui (y, +1, MPFR_RNDN);
186 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
187 if (!mpfr_nan_p (r) || inexact != 0)
188 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN);
189
190 /* fmod (+inf, +1) is NaN */
191 mpfr_neg (x, x, MPFR_RNDN);
192 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
193 if (!mpfr_nan_p (r) || inexact != 0)
194 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN);
195
196 /* fmod (+inf, -inf) is NaN */
197 mpfr_set_inf (y, -1);
198 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
199 if (!mpfr_nan_p (r) || inexact != 0)
200 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN);
201
202 /* fmod (-inf, -inf) is NaN */
203 mpfr_neg (x, x, MPFR_RNDN);
204 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
205 if (!mpfr_nan_p (r) || inexact != 0)
206 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN);
207
208 /* fmod (-inf, +inf) is NaN */
209 mpfr_neg (y, y, MPFR_RNDN);
210 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
211 if (!mpfr_nan_p (r) || inexact != 0)
212 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN);
213
214 /* fmod (+inf, +inf) is NaN */
215 mpfr_neg (x, x, MPFR_RNDN);
216 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
217 if (!mpfr_nan_p (r) || inexact != 0)
218 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN);
219
220 /* fmod (x, +inf) = x, if x is finite */
221 mpfr_set_ui (x, 1, MPFR_RNDN);
222 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
223 if (!mpfr_equal_p (r, x) || inexact != 0)
224 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN);
225 mpfr_neg (x, x, MPFR_RNDN);
226 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
227 if (!mpfr_equal_p (r, x) || inexact != 0)
228 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN);
229
230 /* fmod (+0, +inf) = +0 */
231 mpfr_set_ui (x, 0, MPFR_RNDN);
232 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
233 if (!mpfr_equal_p (r, x) || inexact != 0)
234 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN);
235
236 /* fmod (-0, +inf) = -0 */
237 mpfr_neg (x, x, MPFR_RNDN);
238 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
239 if (!mpfr_equal_p (r, x) || inexact != 0)
240 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN);
241
242 /* fmod (x, -inf) = x, if x is finite */
243 mpfr_set_inf (y, -1);
244 mpfr_set_ui (x, 1, MPFR_RNDN);
245 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
246 if (!mpfr_equal_p (r, x) || inexact != 0)
247 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN);
248 mpfr_neg (x, x, MPFR_RNDN);
249 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
250 if (!mpfr_equal_p (r, x) || inexact != 0)
251 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN);
252
253 /* fmod (+0, -inf) = +0 */
254 mpfr_set_ui (x, 0, MPFR_RNDN);
255 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
256 if (!mpfr_equal_p (r, x) || inexact != 0)
257 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN);
258
259 /* fmod (-0, -inf) = -0 */
260 mpfr_neg (x, x, MPFR_RNDN);
261 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
262 if (!mpfr_equal_p (r, x) || inexact != 0)
263 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN);
264
265 /* fmod (+0, +0) is NaN */
266 mpfr_set_ui (x, 0, MPFR_RNDN);
267 mpfr_set_ui (y, 0, MPFR_RNDN);
268 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
269 if (!mpfr_nan_p (r) || inexact != 0)
270 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN);
271
272 /* fmod (+0, -0) is NaN */
273 mpfr_neg (y, y, MPFR_RNDN);
274 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
275 if (!mpfr_nan_p (r) || inexact != 0)
276 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN);
277
278 /* fmod (+0, +1) = +0 */
279 mpfr_set_ui (y, 1, MPFR_RNDN);
280 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
281 if (!mpfr_equal_p (r, x) || inexact != 0)
282 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN);
283
284 /* fmod (+0, -1) = +0 */
285 mpfr_neg (y, y, MPFR_RNDN);
286 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
287 if (!mpfr_equal_p (r, x) || inexact != 0)
288 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN);
289
290 /* fmod (-0, -1) = -0 */
291 mpfr_neg (x, x, MPFR_RNDN);
292 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
293 if (!mpfr_equal_p (r, x) || inexact != 0)
294 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN);
295
296 /* fmod (-0, +1) = -0 */
297 mpfr_neg (y, y, MPFR_RNDN);
298 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
299 if (!mpfr_equal_p (r, x) || inexact != 0)
300 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN);
301
302 mpfr_set_prec (x, 380);
303 mpfr_set_prec (y, 385);
304 mpfr_set_str_binary (x, "0.11011010010110011101011000100100101100101011010001011100110001100101111001010100001011111110111100101110101010110011010101000100000100011101101100001011101110100111101111111010001001000010000110010110011100111000001110111010000100101001010111100100010001101001110100011110010000000001110001111001101100111011001000110110011100100011111110010100011001000001001011010111010000000000E0");
305 mpfr_set_str_binary (y, "0.1100011000011101011010001100010111001110110111001101010010111100111100011010010011011101111101111001010111111110001001100001111101001000000010100101111001001110010110000111001000101010111001001000100101011111000010100110001111000110011011010101111101100110010101011010011101100001011101001000101111110110110110000001001101110111110110111110111111001001011110001110011111100000000000000E-1");
306 mpfr_set_prec (r, 2);
307 inexact = mpfr_fmod (r, x, y, MPFR_RNDA);
308 mpfr_set_prec (t, 2);
309 mpfr_set_ui_2exp (t, 3, -5, MPFR_RNDN);
310 if (mpfr_cmp_ui_2exp (r, 3, -5) || inexact <= 0)
311 test_failed (r, t, 1, inexact, x, y, MPFR_RNDA);
312
313 mpfr_clears (x, y, r, t, (mpfr_ptr) 0);
314 return;
315 }
316
317 /* bug reported by Eric Veach */
318 static void
bug20090519(void)319 bug20090519 (void)
320 {
321 mpfr_t x, y, r;
322 int inexact;
323
324 mpfr_inits2 (100, x, y, r, (mpfr_ptr) 0);
325
326 mpfr_set_prec (x, 3);
327 mpfr_set_prec (y, 3);
328 mpfr_set_prec (r, 3);
329 mpfr_set_si (x, 8, MPFR_RNDN);
330 mpfr_set_si (y, 7, MPFR_RNDN);
331 check (r, x, y, MPFR_RNDN);
332
333 mpfr_set_prec (x, 10);
334 mpfr_set_prec (y, 10);
335 mpfr_set_prec (r, 10);
336 mpfr_set_ui_2exp (x, 3, 26, MPFR_RNDN);
337 mpfr_set_si (y, (1 << 9) - 1, MPFR_RNDN);
338 check (r, x, y, MPFR_RNDN);
339
340 mpfr_set_prec (x, 100);
341 mpfr_set_prec (y, 100);
342 mpfr_set_prec (r, 100);
343 mpfr_set_str (x, "3.5", 10, MPFR_RNDN);
344 mpfr_set_str (y, "1.1", 10, MPFR_RNDN);
345 check (r, x, y, MPFR_RNDN);
346 /* double check, with a pre-computed value */
347 {
348 mpfr_t er;
349 mpfr_init2 (er, 100);
350 mpfr_set_str (er, "CCCCCCCCCCCCCCCCCCCCCCCC8p-102", 16, MPFR_RNDN);
351
352 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
353 if (!mpfr_equal_p (r, er) || inexact != 0)
354 test_failed (er, r, 0, inexact, x, y, MPFR_RNDN);
355
356 mpfr_clear (er);
357 }
358
359 mpfr_set_si (x, 20, MPFR_RNDN);
360 mpfr_set_ui_2exp (y, 1, 1, MPFR_RNDN); /* exact */
361 mpfr_sin (y, y, MPFR_RNDN);
362 check (r, x, y, MPFR_RNDN);
363
364 mpfr_clears (x, y, r, (mpfr_ptr) 0);
365 }
366
367 static void
bug20160217(void)368 bug20160217 (void)
369 {
370 mpfr_t x, y, r;
371 int inexact, i;
372 mpfr_exp_t emin, emax;
373
374 mpfr_inits2 (53, x, y, r, (mpfr_ptr) 0);
375
376 emin = mpfr_get_emin ();
377 emax = mpfr_get_emax ();
378
379 for (i = 0; i <= 1; i++)
380 {
381 mpfr_set_zero (x, 1);
382 mpfr_nextabove (x);
383 mpfr_set_inf (y, 1);
384 mpfr_nextbelow (y);
385 inexact = mpfr_fmod (r, x, y, MPFR_RNDN);
386 if (!mpfr_equal_p (r, x) || inexact != 0)
387 {
388 printf ("Error for mpfr_fmod (r, nextabove(0), nextbelow(+inf),"
389 " MPFR_RNDN)%s\n", i ? "extended exponent range" : "");
390 printf ("Expected inex = 0, r = ");
391 mpfr_dump (x);
392 printf ("Got inex = %d, r = ", inexact);
393 mpfr_dump (r);
394 exit (1);
395 }
396 set_emin (MPFR_EMIN_MIN);
397 set_emax (MPFR_EMAX_MAX);
398 }
399
400 set_emin (emin);
401 set_emax (emax);
402
403 mpfr_clears (x, y, r, (mpfr_ptr) 0);
404 }
405
406 int
main(int argc,char * argv[])407 main (int argc, char *argv[])
408 {
409 tests_start_mpfr ();
410
411 bug20090519 ();
412 bug20160217 ();
413
414 test_generic (MPFR_PREC_MIN, 100, 100);
415
416 special ();
417 regular ();
418
419 tests_end_mpfr ();
420 return 0;
421 }
422