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