1 /* Test file for mpfr_add_[q,z], mpfr_sub_[q,z], mpfr_div_[q,z],
2 mpfr_mul_[q,z], mpfr_cmp_[f,q,z]
3
4 Copyright 2004-2023 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramba projects, INRIA.
6
7 This file is part of the GNU MPFR Library.
8
9 The GNU MPFR Library is free software; you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
13
14 The GNU MPFR Library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
21 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
22 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23
24 #include "mpfr-test.h"
25
26 #ifndef MPFR_USE_MINI_GMP
27
28 #define CHECK_FOR(str, cond) \
29 if ((cond) == 0) { \
30 printf ("Special case error %s. Ternary value = %d, flags = %u\n", \
31 str, res, __gmpfr_flags); \
32 printf ("Got "); mpfr_dump (y); \
33 printf ("X = "); mpfr_dump (x); \
34 printf ("Q = "); mpz_out_str (stdout, 10, mpq_numref(q)); \
35 printf ("\n /"); mpz_out_str (stdout, 10, mpq_denref(q)); \
36 printf ("\n"); \
37 exit (1); \
38 }
39
40 #define CHECK_FORZ(str, cond) \
41 if ((cond) == 0) { \
42 printf ("Special case error %s. Ternary value = %d, flags = %u\n", \
43 str, res, __gmpfr_flags); \
44 printf ("Got "); mpfr_dump (y); \
45 printf ("X = "); mpfr_dump (x); \
46 printf ("Z = "); mpz_out_str (stdout, 10, z); \
47 printf ("\n"); \
48 exit (1); \
49 }
50
51 static void
special(void)52 special (void)
53 {
54 mpfr_t x, y;
55 mpq_t q;
56 mpz_t z;
57 int res = 0;
58
59 mpfr_init (x);
60 mpfr_init (y);
61 mpq_init (q);
62 mpz_init (z);
63
64 /* cancellation in mpfr_add_q */
65 mpfr_set_prec (x, 60);
66 mpfr_set_prec (y, 20);
67 mpz_set_str (mpq_numref (q), "-187207494", 10);
68 mpz_set_str (mpq_denref (q), "5721", 10);
69 mpfr_set_str_binary (x, "11111111101001011011100101100011011110010011100010000100001E-44");
70 mpfr_add_q (y, x, q, MPFR_RNDN);
71 CHECK_FOR ("cancellation in add_q", mpfr_cmp_ui_2exp (y, 256783, -64) == 0);
72
73 mpfr_set_prec (x, 19);
74 mpfr_set_str_binary (x, "0.1011110101110011100E0");
75 mpz_set_str (mpq_numref (q), "187207494", 10);
76 mpz_set_str (mpq_denref (q), "5721", 10);
77 mpfr_set_prec (y, 29);
78 mpfr_add_q (y, x, q, MPFR_RNDD);
79 mpfr_set_prec (x, 29);
80 mpfr_set_str_binary (x, "11111111101001110011010001001E-14");
81 CHECK_FOR ("cancellation in add_q", mpfr_cmp (x,y) == 0);
82
83 /* Inf */
84 mpfr_set_inf (x, 1);
85 mpz_set_str (mpq_numref (q), "395877315", 10);
86 mpz_set_str (mpq_denref (q), "3508975966", 10);
87 mpfr_set_prec (y, 118);
88 mpfr_add_q (y, x, q, MPFR_RNDU);
89 CHECK_FOR ("inf", mpfr_inf_p (y) && mpfr_sgn (y) > 0);
90 mpfr_sub_q (y, x, q, MPFR_RNDU);
91 CHECK_FOR ("inf", mpfr_inf_p (y) && mpfr_sgn (y) > 0);
92
93 /* Nan */
94 MPFR_SET_NAN (x);
95 mpfr_add_q (y, x, q, MPFR_RNDU);
96 CHECK_FOR ("nan", mpfr_nan_p (y));
97 mpfr_sub_q (y, x, q, MPFR_RNDU);
98 CHECK_FOR ("nan", mpfr_nan_p (y));
99
100 /* Exact value */
101 mpfr_set_prec (x, 60);
102 mpfr_set_prec (y, 60);
103 mpfr_set_str1 (x, "0.5");
104 mpz_set_str (mpq_numref (q), "3", 10);
105 mpz_set_str (mpq_denref (q), "2", 10);
106 res = mpfr_add_q (y, x, q, MPFR_RNDU);
107 CHECK_FOR ("0.5+3/2", mpfr_cmp_ui(y, 2)==0 && res==0);
108 res = mpfr_sub_q (y, x, q, MPFR_RNDU);
109 CHECK_FOR ("0.5-3/2", mpfr_cmp_si(y, -1)==0 && res==0);
110
111 /* Inf rational */
112 mpq_set_ui (q, 1, 0);
113 mpfr_set_str1 (x, "0.5");
114 res = mpfr_add_q (y, x, q, MPFR_RNDN);
115 CHECK_FOR ("0.5+1/0", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0);
116 res = mpfr_sub_q (y, x, q, MPFR_RNDN);
117 CHECK_FOR ("0.5-1/0", mpfr_inf_p (y) && MPFR_IS_NEG (y) && res == 0);
118 mpq_set_si (q, -1, 0);
119 res = mpfr_add_q (y, x, q, MPFR_RNDN);
120 CHECK_FOR ("0.5+ -1/0", mpfr_inf_p (y) && MPFR_IS_NEG (y) && res == 0);
121 res = mpfr_sub_q (y, x, q, MPFR_RNDN);
122 CHECK_FOR ("0.5- -1/0", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0);
123 res = mpfr_div_q (y, x, q, MPFR_RNDN);
124 CHECK_FOR ("0.5 / (-1/0)", mpfr_zero_p (y) && MPFR_IS_NEG (y) && res == 0);
125 mpq_set_ui (q, 1, 0);
126 mpfr_set_inf (x, 1);
127 res = mpfr_add_q (y, x, q, MPFR_RNDN);
128 CHECK_FOR ("+Inf + +Inf", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0);
129 res = mpfr_sub_q (y, x, q, MPFR_RNDN);
130 CHECK_FOR ("+Inf - +Inf", MPFR_IS_NAN (y) && res == 0);
131 mpfr_set_inf (x, -1);
132 res = mpfr_add_q (y, x, q, MPFR_RNDN);
133 CHECK_FOR ("-Inf + +Inf", MPFR_IS_NAN (y) && res == 0);
134 res = mpfr_sub_q (y, x, q, MPFR_RNDN);
135 CHECK_FOR ("-Inf - +Inf", mpfr_inf_p (y) && MPFR_IS_NEG (y) && res == 0);
136 mpq_set_si (q, -1, 0);
137 mpfr_set_inf (x, 1);
138 res = mpfr_add_q (y, x, q, MPFR_RNDN);
139 CHECK_FOR ("+Inf + -Inf", MPFR_IS_NAN (y) && res == 0);
140 res = mpfr_sub_q (y, x, q, MPFR_RNDN);
141 CHECK_FOR ("+Inf - -Inf", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0);
142 mpfr_set_inf (x, -1);
143 res = mpfr_add_q (y, x, q, MPFR_RNDN);
144 CHECK_FOR ("-Inf + -Inf", mpfr_inf_p (y) && MPFR_IS_NEG (y) && res == 0);
145 res = mpfr_sub_q (y, x, q, MPFR_RNDN);
146 CHECK_FOR ("-Inf - -Inf", MPFR_IS_NAN (y) && res == 0);
147
148 /* 0 */
149 mpq_set_ui (q, 0, 1);
150 mpfr_set_ui (x, 42, MPFR_RNDN);
151 res = mpfr_add_q (y, x, q, MPFR_RNDN);
152 CHECK_FOR ("42+0/1", mpfr_cmp_ui (y, 42) == 0 && res == 0);
153 res = mpfr_sub_q (y, x, q, MPFR_RNDN);
154 CHECK_FOR ("42-0/1", mpfr_cmp_ui (y, 42) == 0 && res == 0);
155 res = mpfr_mul_q (y, x, q, MPFR_RNDN);
156 CHECK_FOR ("42*0/1", mpfr_zero_p (y) && MPFR_IS_POS (y) && res == 0);
157 mpfr_clear_flags ();
158 res = mpfr_div_q (y, x, q, MPFR_RNDN);
159 CHECK_FOR ("42/(0/1)", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0
160 && mpfr_divby0_p ());
161 mpz_set_ui (z, 0);
162 mpfr_clear_flags ();
163 res = mpfr_div_z (y, x, z, MPFR_RNDN);
164 CHECK_FORZ ("42/0", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0
165 && mpfr_divby0_p ());
166
167 mpz_clear (z);
168 mpq_clear (q);
169 mpfr_clear (x);
170 mpfr_clear (y);
171 }
172
173 static void
check_for_zero(void)174 check_for_zero (void)
175 {
176 /* Check that 0 is unsigned! */
177 mpq_t q;
178 mpz_t z;
179 mpfr_t x;
180 int r;
181 mpfr_sign_t i;
182
183 mpfr_init (x);
184 mpz_init (z);
185 mpq_init (q);
186
187 mpz_set_ui (z, 0);
188 mpq_set_ui (q, 0, 1);
189
190 MPFR_SET_ZERO (x);
191 RND_LOOP (r)
192 {
193 for (i = MPFR_SIGN_NEG ; i <= MPFR_SIGN_POS ;
194 i+=MPFR_SIGN_POS-MPFR_SIGN_NEG)
195 {
196 MPFR_SET_SIGN(x, i);
197 mpfr_add_z (x, x, z, (mpfr_rnd_t) r);
198 if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
199 {
200 printf("GMP Zero errors for add_z & rnd=%s & s=%d\n",
201 mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
202 mpfr_dump (x);
203 exit (1);
204 }
205 mpfr_sub_z (x, x, z, (mpfr_rnd_t) r);
206 if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
207 {
208 printf("GMP Zero errors for sub_z & rnd=%s & s=%d\n",
209 mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
210 mpfr_dump (x);
211 exit (1);
212 }
213 mpfr_mul_z (x, x, z, (mpfr_rnd_t) r);
214 if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
215 {
216 printf("GMP Zero errors for mul_z & rnd=%s & s=%d\n",
217 mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
218 mpfr_dump (x);
219 exit (1);
220 }
221 mpfr_add_q (x, x, q, (mpfr_rnd_t) r);
222 if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
223 {
224 printf("GMP Zero errors for add_q & rnd=%s & s=%d\n",
225 mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
226 mpfr_dump (x);
227 exit (1);
228 }
229 mpfr_sub_q (x, x, q, (mpfr_rnd_t) r);
230 if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
231 {
232 printf("GMP Zero errors for sub_q & rnd=%s & s=%d\n",
233 mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
234 mpfr_dump (x);
235 exit (1);
236 }
237 }
238 }
239
240 mpq_clear (q);
241 mpz_clear (z);
242 mpfr_clear (x);
243 }
244
245 static void
test_cmp_z(mpfr_prec_t pmin,mpfr_prec_t pmax,int nmax)246 test_cmp_z (mpfr_prec_t pmin, mpfr_prec_t pmax, int nmax)
247 {
248 mpfr_t x, z;
249 mpz_t y;
250 mpfr_prec_t p;
251 int res1, res2;
252 int n;
253
254 mpfr_init (x);
255 mpfr_init2 (z, MPFR_PREC_MIN);
256 mpz_init (y);
257
258 /* check the erange flag when x is NaN */
259 mpfr_set_nan (x);
260 mpz_set_ui (y, 17);
261 mpfr_clear_erangeflag ();
262 res1 = mpfr_cmp_z (x, y);
263 if (res1 != 0 || mpfr_erangeflag_p () == 0)
264 {
265 printf ("Error for mpfr_cmp_z (NaN, 17)\n");
266 printf ("Return value: expected 0, got %d\n", res1);
267 printf ("Erange flag: expected set, got %d\n", mpfr_erangeflag_p ());
268 exit (1);
269 }
270
271 for (p = pmin ; p < pmax ; p++)
272 {
273 mpfr_set_prec (x, p);
274 for (n = 0 ; n < nmax ; n++)
275 {
276 mpfr_urandomb (x, RANDS);
277 mpz_urandomb (y, RANDS, 1024);
278 if (!MPFR_IS_SINGULAR (x))
279 {
280 mpfr_sub_z (z, x, y, MPFR_RNDN);
281 res1 = (mpfr_sgn) (z);
282 res2 = mpfr_cmp_z (x, y);
283 if (res1 != res2)
284 {
285 printf("Error for mpfr_cmp_z: res=%d sub_z gives %d\n",
286 res2, res1);
287 exit (1);
288 }
289 }
290 }
291 }
292 mpz_clear (y);
293 mpfr_clear (x);
294 mpfr_clear (z);
295 }
296
297 static void
test_cmp_q(mpfr_prec_t pmin,mpfr_prec_t pmax,int nmax)298 test_cmp_q (mpfr_prec_t pmin, mpfr_prec_t pmax, int nmax)
299 {
300 mpfr_t x, z;
301 mpq_t y;
302 mpfr_prec_t p;
303 int res1, res2;
304 int n;
305
306 mpfr_init (x);
307 mpfr_init2 (z, MPFR_PREC_MIN);
308 mpq_init (y);
309
310 /* Check the flags when x is NaN: the erange flags must be set, and
311 only this one. */
312 mpfr_set_nan (x);
313 mpq_set_ui (y, 17, 1);
314 mpfr_clear_flags ();
315 res1 = mpfr_cmp_q (x, y);
316 if (res1 != 0 || __gmpfr_flags != MPFR_FLAGS_ERANGE)
317 {
318 printf ("Error for mpfr_cmp_q (NaN, 17)\n");
319 printf ("Return value: expected 0, got %d\n", res1);
320 printf ("Expected flags:");
321 flags_out (MPFR_FLAGS_ERANGE);
322 printf ("Got flags: ");
323 flags_out (__gmpfr_flags);
324 exit (1);
325 }
326
327 /* Check the flags when y is NaN: the erange flags must be set, and
328 only this one. */
329 mpfr_set_ui (x, 42, MPFR_RNDN);
330 /* A NaN rational is represented by 0/0 (MPFR extension). */
331 mpz_set_ui (mpq_numref (y), 0);
332 mpz_set_ui (mpq_denref (y), 0);
333 mpfr_clear_flags ();
334 res1 = mpfr_cmp_q (x, y);
335 if (res1 != 0 || __gmpfr_flags != MPFR_FLAGS_ERANGE)
336 {
337 printf ("Error for mpfr_cmp_q (42, NaN)\n");
338 printf ("Return value: expected 0, got %d\n", res1);
339 printf ("Expected flags:");
340 flags_out (MPFR_FLAGS_ERANGE);
341 printf ("Got flags: ");
342 flags_out (__gmpfr_flags);
343 exit (1);
344 }
345
346 for (p = pmin ; p < pmax ; p++)
347 {
348 mpfr_set_prec (x, p);
349 for (n = 0 ; n < nmax ; n++)
350 {
351 mpfr_urandomb (x, RANDS);
352 mpq_set_ui (y, randlimb (), randlimb() );
353 if (!MPFR_IS_SINGULAR (x))
354 {
355 mpfr_sub_q (z, x, y, MPFR_RNDN);
356 res1 = (mpfr_sgn) (z);
357 res2 = mpfr_cmp_q (x, y);
358 if (res1 != res2)
359 {
360 printf("Error for mpfr_cmp_q: res=%d sub_z gives %d\n",
361 res2, res1);
362 exit (1);
363 }
364 }
365 }
366 }
367
368 /* check for y = 1/0 */
369 mpz_set_ui (mpq_numref (y), 1);
370 mpz_set_ui (mpq_denref (y), 0);
371 mpfr_set_ui (x, 1, MPFR_RNDN);
372 MPFR_ASSERTN(mpfr_cmp_q (x, y) < 0);
373 mpfr_set_inf (x, -1);
374 MPFR_ASSERTN(mpfr_cmp_q (x, y) < 0);
375 mpfr_set_inf (x, +1);
376 MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0);
377 mpfr_set_nan (x);
378 mpfr_clear_erangeflag ();
379 MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0);
380 MPFR_ASSERTN(mpfr_erangeflag_p ());
381
382 /* check for y = -1/0 */
383 mpz_set_si (mpq_numref (y), -1);
384 mpz_set_ui (mpq_denref (y), 0);
385 mpfr_set_ui (x, 1, MPFR_RNDN);
386 MPFR_ASSERTN(mpfr_cmp_q (x, y) > 0);
387 mpfr_set_inf (x, -1);
388 MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0);
389 mpfr_set_inf (x, +1);
390 MPFR_ASSERTN(mpfr_cmp_q (x, y) > 0);
391 mpfr_set_nan (x);
392 mpfr_clear_erangeflag ();
393 MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0);
394 MPFR_ASSERTN(mpfr_erangeflag_p ());
395
396 /* check for y = 0/0 */
397 mpz_set_ui (mpq_numref (y), 0);
398 mpz_set_ui (mpq_denref (y), 0);
399 mpfr_set_ui (x, 1, MPFR_RNDN);
400 mpfr_clear_erangeflag ();
401 MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0);
402 MPFR_ASSERTN(mpfr_erangeflag_p ());
403 mpfr_set_inf (x, -1);
404 mpfr_clear_erangeflag ();
405 MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0);
406 MPFR_ASSERTN(mpfr_erangeflag_p ());
407 mpfr_set_inf (x, +1);
408 mpfr_clear_erangeflag ();
409 MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0);
410 MPFR_ASSERTN(mpfr_erangeflag_p ());
411 mpfr_set_nan (x);
412 mpfr_clear_erangeflag ();
413 MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0);
414 MPFR_ASSERTN(mpfr_erangeflag_p ());
415
416 mpq_clear (y);
417 mpfr_clear (x);
418 mpfr_clear (z);
419 }
420
421 static void
test_cmp_f(mpfr_prec_t pmin,mpfr_prec_t pmax,int nmax)422 test_cmp_f (mpfr_prec_t pmin, mpfr_prec_t pmax, int nmax)
423 {
424 mpfr_t x, z;
425 mpf_t y;
426 mpfr_prec_t p;
427 int res1, res2;
428 int n;
429
430 mpfr_init (x);
431 mpfr_init2 (z, pmax+GMP_NUMB_BITS);
432 mpf_init2 (y, MPFR_PREC_MIN);
433
434 /* check the erange flag when x is NaN */
435 mpfr_set_nan (x);
436 mpf_set_ui (y, 17);
437 mpfr_clear_erangeflag ();
438 res1 = mpfr_cmp_f (x, y);
439 if (res1 != 0 || mpfr_erangeflag_p () == 0)
440 {
441 printf ("Error for mpfr_cmp_f (NaN, 17)\n");
442 printf ("Return value: expected 0, got %d\n", res1);
443 printf ("Erange flag: expected set, got %d\n", mpfr_erangeflag_p ());
444 exit (1);
445 }
446
447 for (p = pmin ; p < pmax ; p += 3)
448 {
449 mpfr_set_prec (x, p);
450 mpf_set_prec (y, p);
451 for (n = 0 ; n < nmax ; n++)
452 {
453 mpfr_urandomb (x, RANDS);
454 mpf_urandomb (y, RANDS, p);
455 if (!MPFR_IS_SINGULAR (x))
456 {
457 mpfr_set_f (z, y, MPFR_RNDN);
458 mpfr_sub (z, x, z, MPFR_RNDN);
459 res1 = (mpfr_sgn) (z);
460 res2 = mpfr_cmp_f (x, y);
461 if (res1 != res2)
462 {
463 printf("Error for mpfr_cmp_f: res=%d sub gives %d\n",
464 res2, res1);
465 exit (1);
466 }
467 }
468 }
469 }
470 mpf_clear (y);
471 mpfr_clear (x);
472 mpfr_clear (z);
473 }
474
475 static void
test_specialz(int (* mpfr_func)(mpfr_ptr,mpfr_srcptr,mpz_srcptr,mpfr_rnd_t),void (* mpz_func)(mpz_ptr,mpz_srcptr,mpz_srcptr),const char * op)476 test_specialz (int (*mpfr_func)(mpfr_ptr, mpfr_srcptr, mpz_srcptr, mpfr_rnd_t),
477 void (*mpz_func)(mpz_ptr, mpz_srcptr, mpz_srcptr),
478 const char *op)
479 {
480 mpfr_t x1, x2;
481 mpz_t z1, z2;
482 int res;
483
484 mpfr_inits2 (128, x1, x2, (mpfr_ptr) 0);
485 mpz_init (z1); mpz_init(z2);
486 mpz_fac_ui (z1, 19); /* 19!+1 fits perfectly in a 128 bits mantissa */
487 mpz_add_ui (z1, z1, 1);
488 mpz_fac_ui (z2, 20); /* 20!+1 fits perfectly in a 128 bits mantissa */
489 mpz_add_ui (z2, z2, 1);
490
491 res = mpfr_set_z(x1, z1, MPFR_RNDN);
492 if (res)
493 {
494 printf("Specialz %s: set_z1 error\n", op);
495 exit(1);
496 }
497 mpfr_set_z (x2, z2, MPFR_RNDN);
498 if (res)
499 {
500 printf("Specialz %s: set_z2 error\n", op);
501 exit(1);
502 }
503
504 /* (19!+1) * (20!+1) fits in a 128 bits number */
505 res = mpfr_func(x1, x1, z2, MPFR_RNDN);
506 if (res)
507 {
508 printf("Specialz %s: wrong inexact flag.\n", op);
509 exit(1);
510 }
511 mpz_func(z1, z1, z2);
512 res = mpfr_set_z (x2, z1, MPFR_RNDN);
513 if (res)
514 {
515 printf("Specialz %s: set_z2 error\n", op);
516 exit(1);
517 }
518 if (mpfr_cmp(x1, x2))
519 {
520 printf("Specialz %s: results differ.\nx1=", op);
521 mpfr_dump (x1);
522 printf ("x2=");
523 mpfr_dump (x2);
524 printf ("Z2=");
525 mpz_out_str (stdout, 2, z1);
526 putchar('\n');
527 exit(1);
528 }
529
530 mpz_set_ui (z1, 1);
531 mpz_set_ui (z2, 0);
532 mpfr_set_ui (x1, 1, MPFR_RNDN);
533 mpz_func (z1, z1, z2);
534 res = mpfr_func(x1, x1, z2, MPFR_RNDN);
535 mpfr_set_z (x2, z1, MPFR_RNDN);
536 if (mpfr_cmp(x1, x2))
537 {
538 printf("Specialz %s: results differ(2).\nx1=", op);
539 mpfr_dump (x1);
540 printf ("x2=");
541 mpfr_dump (x2);
542 exit(1);
543 }
544
545 mpz_clear (z1); mpz_clear(z2);
546 mpfr_clears (x1, x2, (mpfr_ptr) 0);
547 }
548
549 static void
test_special2z(int (* mpfr_func)(mpfr_ptr,mpz_srcptr,mpfr_srcptr,mpfr_rnd_t),void (* mpz_func)(mpz_ptr,mpz_srcptr,mpz_srcptr),const char * op)550 test_special2z (int (*mpfr_func)(mpfr_ptr, mpz_srcptr, mpfr_srcptr, mpfr_rnd_t),
551 void (*mpz_func)(mpz_ptr, mpz_srcptr, mpz_srcptr),
552 const char *op)
553 {
554 mpfr_t x1, x2;
555 mpz_t z1, z2;
556 int res;
557
558 mpfr_inits2 (128, x1, x2, (mpfr_ptr) 0);
559 mpz_init (z1); mpz_init(z2);
560 mpz_fac_ui (z1, 19); /* 19!+1 fits perfectly in a 128 bits mantissa */
561 mpz_add_ui (z1, z1, 1);
562 mpz_fac_ui (z2, 20); /* 20!+1 fits perfectly in a 128 bits mantissa */
563 mpz_add_ui (z2, z2, 1);
564
565 res = mpfr_set_z(x1, z1, MPFR_RNDN);
566 if (res)
567 {
568 printf("Special2z %s: set_z1 error\n", op);
569 exit(1);
570 }
571 mpfr_set_z (x2, z2, MPFR_RNDN);
572 if (res)
573 {
574 printf("Special2z %s: set_z2 error\n", op);
575 exit(1);
576 }
577
578 /* (19!+1) * (20!+1) fits in a 128 bits number */
579 res = mpfr_func(x1, z1, x2, MPFR_RNDN);
580 if (res)
581 {
582 printf("Special2z %s: wrong inexact flag.\n", op);
583 exit(1);
584 }
585 mpz_func(z1, z1, z2);
586 res = mpfr_set_z (x2, z1, MPFR_RNDN);
587 if (res)
588 {
589 printf("Special2z %s: set_z2 error\n", op);
590 exit(1);
591 }
592 if (mpfr_cmp(x1, x2))
593 {
594 printf("Special2z %s: results differ.\nx1=", op);
595 mpfr_dump (x1);
596 printf ("x2=");
597 mpfr_dump (x2);
598 printf ("Z2=");
599 mpz_out_str (stdout, 2, z1);
600 putchar('\n');
601 exit(1);
602 }
603
604 mpz_set_ui (z1, 0);
605 mpz_set_ui (z2, 1);
606 mpfr_set_ui (x2, 1, MPFR_RNDN);
607 res = mpfr_func(x1, z1, x2, MPFR_RNDN);
608 mpz_func (z1, z1, z2);
609 mpfr_set_z (x2, z1, MPFR_RNDN);
610 if (mpfr_cmp(x1, x2))
611 {
612 printf("Special2z %s: results differ(2).\nx1=", op);
613 mpfr_dump (x1);
614 printf ("x2=");
615 mpfr_dump (x2);
616 exit(1);
617 }
618
619 mpz_clear (z1); mpz_clear(z2);
620 mpfr_clears (x1, x2, (mpfr_ptr) 0);
621 }
622
623 static void
test_genericz(mpfr_prec_t p0,mpfr_prec_t p1,unsigned int N,int (* func)(mpfr_ptr,mpfr_srcptr,mpz_srcptr,mpfr_rnd_t),const char * op)624 test_genericz (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int N,
625 int (*func)(mpfr_ptr, mpfr_srcptr, mpz_srcptr, mpfr_rnd_t),
626 const char *op)
627 {
628 mpfr_prec_t prec;
629 mpfr_t arg1, dst_big, dst_small, tmp;
630 mpz_t arg2;
631 mpfr_rnd_t rnd;
632 int inexact, compare, compare2;
633 unsigned int n;
634
635 mpfr_inits (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0);
636 mpz_init (arg2);
637
638 for (prec = p0; prec <= p1; prec++)
639 {
640 mpfr_set_prec (arg1, prec);
641 mpfr_set_prec (tmp, prec);
642 mpfr_set_prec (dst_small, prec);
643
644 for (n=0; n<N; n++)
645 {
646 mpfr_urandomb (arg1, RANDS);
647 mpz_urandomb (arg2, RANDS, 1024);
648 rnd = RND_RAND_NO_RNDF ();
649 mpfr_set_prec (dst_big, 2*prec);
650 compare = func (dst_big, arg1, arg2, rnd);
651 if (mpfr_can_round (dst_big, 2*prec, rnd, rnd, prec))
652 {
653 mpfr_set (tmp, dst_big, rnd);
654 inexact = func (dst_small, arg1, arg2, rnd);
655 if (mpfr_cmp (tmp, dst_small))
656 {
657 printf ("Results differ for prec=%u rnd_mode=%s and %s_z:\n"
658 "arg1=",
659 (unsigned) prec, mpfr_print_rnd_mode (rnd), op);
660 mpfr_dump (arg1);
661 printf ("arg2=");
662 mpz_out_str (stdout, 10, arg2);
663 printf ("\ngot ");
664 mpfr_dump (dst_small);
665 printf ("expected ");
666 mpfr_dump (tmp);
667 printf ("approx ");
668 mpfr_dump (dst_big);
669 exit (1);
670 }
671 compare2 = mpfr_cmp (tmp, dst_big);
672 /* if rounding to nearest, cannot know the sign of t - f(x)
673 because of composed rounding: y = o(f(x)) and t = o(y) */
674 if (compare * compare2 >= 0)
675 compare = compare + compare2;
676 else
677 compare = inexact; /* cannot determine sign(t-f(x)) */
678 if (((inexact == 0) && (compare != 0)) ||
679 ((inexact > 0) && (compare <= 0)) ||
680 ((inexact < 0) && (compare >= 0)))
681 {
682 printf ("Wrong inexact flag for rnd=%s and %s_z:\n"
683 "expected %d, got %d\n",
684 mpfr_print_rnd_mode (rnd), op, compare, inexact);
685 printf ("arg1="); mpfr_dump (arg1);
686 printf ("arg2="); mpz_out_str(stdout, 2, arg2);
687 printf ("\ndstl="); mpfr_dump (dst_big);
688 printf ("dsts="); mpfr_dump (dst_small);
689 printf ("tmp ="); mpfr_dump (tmp);
690 exit (1);
691 }
692 }
693 }
694 }
695
696 mpz_clear (arg2);
697 mpfr_clears (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0);
698 }
699
700 static void
test_generic2z(mpfr_prec_t p0,mpfr_prec_t p1,unsigned int N,int (* func)(mpfr_ptr,mpz_srcptr,mpfr_srcptr,mpfr_rnd_t),const char * op)701 test_generic2z (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int N,
702 int (*func)(mpfr_ptr, mpz_srcptr, mpfr_srcptr, mpfr_rnd_t),
703 const char *op)
704 {
705 mpfr_prec_t prec;
706 mpfr_t arg1, dst_big, dst_small, tmp;
707 mpz_t arg2;
708 mpfr_rnd_t rnd;
709 int inexact, compare, compare2;
710 unsigned int n;
711
712 mpfr_inits (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0);
713 mpz_init (arg2);
714
715 for (prec = p0; prec <= p1; prec++)
716 {
717 mpfr_set_prec (arg1, prec);
718 mpfr_set_prec (tmp, prec);
719 mpfr_set_prec (dst_small, prec);
720
721 for (n=0; n<N; n++)
722 {
723 mpfr_urandomb (arg1, RANDS);
724 mpz_urandomb (arg2, RANDS, 1024);
725 rnd = RND_RAND_NO_RNDF ();
726 mpfr_set_prec (dst_big, 2*prec);
727 compare = func(dst_big, arg2, arg1, rnd);
728 if (mpfr_can_round (dst_big, 2*prec, rnd, rnd, prec))
729 {
730 mpfr_set (tmp, dst_big, rnd);
731 inexact = func(dst_small, arg2, arg1, rnd);
732 if (mpfr_cmp (tmp, dst_small))
733 {
734 printf ("Results differ for prec=%u rnd_mode=%s and %s_z:\n"
735 "arg1=",
736 (unsigned) prec, mpfr_print_rnd_mode (rnd), op);
737 mpfr_dump (arg1);
738 printf ("arg2=");
739 mpz_out_str (stdout, 10, arg2);
740 printf ("\ngot ");
741 mpfr_dump (dst_small);
742 printf ("expected ");
743 mpfr_dump (tmp);
744 printf ("approx ");
745 mpfr_dump (dst_big);
746 exit (1);
747 }
748 compare2 = mpfr_cmp (tmp, dst_big);
749 /* if rounding to nearest, cannot know the sign of t - f(x)
750 because of composed rounding: y = o(f(x)) and t = o(y) */
751 if (compare * compare2 >= 0)
752 compare = compare + compare2;
753 else
754 compare = inexact; /* cannot determine sign(t-f(x)) */
755 if (((inexact == 0) && (compare != 0)) ||
756 ((inexact > 0) && (compare <= 0)) ||
757 ((inexact < 0) && (compare >= 0)))
758 {
759 printf ("Wrong inexact flag for rnd=%s and %s_z:\n"
760 "expected %d, got %d\n",
761 mpfr_print_rnd_mode (rnd), op, compare, inexact);
762 printf ("arg1="); mpfr_dump (arg1);
763 printf ("arg2="); mpz_out_str(stdout, 2, arg2);
764 printf ("\ndstl="); mpfr_dump (dst_big);
765 printf ("dsts="); mpfr_dump (dst_small);
766 printf ("tmp ="); mpfr_dump (tmp);
767 exit (1);
768 }
769 }
770 }
771 }
772
773 mpz_clear (arg2);
774 mpfr_clears (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0);
775 }
776
777 static void
test_genericq(mpfr_prec_t p0,mpfr_prec_t p1,unsigned int N,int (* func)(mpfr_ptr,mpfr_srcptr,mpq_srcptr,mpfr_rnd_t),const char * op)778 test_genericq (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int N,
779 int (*func)(mpfr_ptr, mpfr_srcptr, mpq_srcptr, mpfr_rnd_t),
780 const char *op)
781 {
782 mpfr_prec_t prec;
783 mpfr_t arg1, dst_big, dst_small, tmp;
784 mpq_t arg2;
785 mpfr_rnd_t rnd;
786 int inexact, compare, compare2;
787 unsigned int n;
788
789 mpfr_inits (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0);
790 mpq_init (arg2);
791
792 for (prec = p0; prec <= p1; prec++)
793 {
794 mpfr_set_prec (arg1, prec);
795 mpfr_set_prec (tmp, prec);
796 mpfr_set_prec (dst_small, prec);
797
798 for (n=0; n<N; n++)
799 {
800 mpfr_urandomb (arg1, RANDS);
801 mpq_set_ui (arg2, randlimb (), randlimb() );
802 mpq_canonicalize (arg2);
803 rnd = RND_RAND_NO_RNDF ();
804 mpfr_set_prec (dst_big, prec+10);
805 compare = func(dst_big, arg1, arg2, rnd);
806 if (mpfr_can_round (dst_big, prec+10, rnd, rnd, prec))
807 {
808 mpfr_set (tmp, dst_big, rnd);
809 inexact = func(dst_small, arg1, arg2, rnd);
810 if (mpfr_cmp (tmp, dst_small))
811 {
812 printf ("Results differ for prec=%u rnd_mode=%s and %s_q:\n"
813 "arg1=",
814 (unsigned) prec, mpfr_print_rnd_mode (rnd), op);
815 mpfr_dump (arg1);
816 printf ("arg2=");
817 mpq_out_str(stdout, 2, arg2);
818 printf ("\ngot ");
819 mpfr_dump (dst_small);
820 printf ("expected ");
821 mpfr_dump (tmp);
822 printf ("approx ");
823 mpfr_dump (dst_big);
824 exit (1);
825 }
826 compare2 = mpfr_cmp (tmp, dst_big);
827 /* if rounding to nearest, cannot know the sign of t - f(x)
828 because of composed rounding: y = o(f(x)) and t = o(y) */
829 if (compare * compare2 >= 0)
830 compare = compare + compare2;
831 else
832 compare = inexact; /* cannot determine sign(t-f(x)) */
833 if (((inexact == 0) && (compare != 0)) ||
834 ((inexact > 0) && (compare <= 0)) ||
835 ((inexact < 0) && (compare >= 0)))
836 {
837 printf ("Wrong inexact flag for rnd=%s and %s_q:\n"
838 "expected %d, got %d",
839 mpfr_print_rnd_mode (rnd), op, compare, inexact);
840 printf ("arg1="); mpfr_dump (arg1);
841 printf ("arg2="); mpq_out_str(stdout, 2, arg2);
842 printf ("\ndstl="); mpfr_dump (dst_big);
843 printf ("dsts="); mpfr_dump (dst_small);
844 printf ("tmp ="); mpfr_dump (tmp);
845 exit (1);
846 }
847 }
848 }
849 }
850
851 mpq_clear (arg2);
852 mpfr_clears (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0);
853 }
854
855 static void
test_specialq(mpfr_prec_t p0,mpfr_prec_t p1,unsigned int N,int (* mpfr_func)(mpfr_ptr,mpfr_srcptr,mpq_srcptr,mpfr_rnd_t),void (* mpq_func)(mpq_ptr,mpq_srcptr,mpq_srcptr),const char * op)856 test_specialq (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int N,
857 int (*mpfr_func)(mpfr_ptr, mpfr_srcptr, mpq_srcptr, mpfr_rnd_t),
858 void (*mpq_func)(mpq_ptr, mpq_srcptr, mpq_srcptr),
859 const char *op)
860 {
861 mpfr_t fra, frb, frq;
862 mpq_t q1, q2, qr;
863 unsigned int n;
864 mpfr_prec_t prec;
865
866 for (prec = p0 ; prec < p1 ; prec++)
867 {
868 mpfr_inits2 (prec, fra, frb, frq, (mpfr_ptr) 0);
869 mpq_init (q1); mpq_init(q2); mpq_init (qr);
870
871 for (n = 0 ; n < N ; n++)
872 {
873 mpq_set_ui(q1, randlimb(), randlimb() );
874 mpq_set_ui(q2, randlimb(), randlimb() );
875 mpq_canonicalize (q1);
876 mpq_canonicalize (q2);
877 mpq_func (qr, q1, q2);
878 mpfr_set_q (fra, q1, MPFR_RNDD);
879 mpfr_func (fra, fra, q2, MPFR_RNDD);
880 mpfr_set_q (frb, q1, MPFR_RNDU);
881 mpfr_func (frb, frb, q2, MPFR_RNDU);
882 mpfr_set_q (frq, qr, MPFR_RNDN);
883 /* We should have fra <= qr <= frb */
884 if ( (mpfr_cmp(fra, frq) > 0) || (mpfr_cmp (frq, frb) > 0))
885 {
886 printf("Range error for prec=%lu and %s",
887 (unsigned long) prec, op);
888 printf ("\nq1="); mpq_out_str(stdout, 2, q1);
889 printf ("\nq2="); mpq_out_str(stdout, 2, q2);
890 printf ("\nfr_dn="); mpfr_dump (fra);
891 printf ("fr_q ="); mpfr_dump (frq);
892 printf ("fr_up="); mpfr_dump (frb);
893 exit (1);
894 }
895 }
896
897 mpq_clear (q1); mpq_clear (q2); mpq_clear (qr);
898 mpfr_clears (fra, frb, frq, (mpfr_ptr) 0);
899 }
900 }
901
902 static void
bug_mul_q_20100810(void)903 bug_mul_q_20100810 (void)
904 {
905 mpfr_t x;
906 mpfr_t y;
907 mpq_t q;
908 int inexact;
909
910 mpfr_init (x);
911 mpfr_init (y);
912 mpq_init (q);
913
914 /* mpfr_mul_q: the inexact value must be set in case of overflow */
915 mpq_set_ui (q, 4096, 3);
916 mpfr_set_inf (x, +1);
917 mpfr_nextbelow (x);
918 inexact = mpfr_mul_q (y, x, q, MPFR_RNDU);
919
920 if (inexact <= 0)
921 {
922 printf ("Overflow error in mpfr_mul_q. ");
923 printf ("Wrong inexact flag: got %d, should be positive.\n", inexact);
924
925 exit (1);
926 }
927 if (!mpfr_inf_p (y))
928 {
929 printf ("Overflow error in mpfr_mul_q (y, x, q, MPFR_RNDD). ");
930 printf ("\nx = ");
931 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDD);
932 printf ("\nq = ");
933 mpq_out_str (stdout, 10, q);
934 printf ("\ny = ");
935 mpfr_out_str (stdout, 10, 0, y, MPFR_RNDD);
936 printf (" (should be +infinity)\n");
937
938 exit (1);
939 }
940
941 mpq_clear (q);
942 mpfr_clear (y);
943 mpfr_clear (x);
944 }
945
946 static void
bug_div_q_20100810(void)947 bug_div_q_20100810 (void)
948 {
949 mpfr_t x;
950 mpfr_t y;
951 mpq_t q;
952 int inexact;
953
954 mpfr_init (x);
955 mpfr_init (y);
956 mpq_init (q);
957
958 /* mpfr_div_q: the inexact value must be set in case of overflow */
959 mpq_set_ui (q, 3, 4096);
960 mpfr_set_inf (x, +1);
961 mpfr_nextbelow (x);
962 inexact = mpfr_div_q (y, x, q, MPFR_RNDU);
963
964 if (inexact <= 0)
965 {
966 printf ("Overflow error in mpfr_div_q. ");
967 printf ("Wrong inexact flag: got %d, should be positive.\n", inexact);
968
969 exit (1);
970 }
971 if (!mpfr_inf_p (y))
972 {
973 printf ("Overflow error in mpfr_div_q (y, x, q, MPFR_RNDD). ");
974 printf ("\nx = ");
975 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDD);
976 printf ("\nq = ");
977 mpq_out_str (stdout, 10, q);
978 printf ("\ny = ");
979 mpfr_out_str (stdout, 10, 0, y, MPFR_RNDD);
980 printf (" (should be +infinity)\n");
981
982 exit (1);
983 }
984
985 mpq_clear (q);
986 mpfr_clear (y);
987 mpfr_clear (x);
988 }
989
990 static void
bug_mul_div_q_20100818(void)991 bug_mul_div_q_20100818 (void)
992 {
993 mpq_t qa, qb;
994 mpfr_t x1, x2, y1, y2, y3;
995 mpfr_exp_t emin, emax, e;
996 int inex;
997 int rnd;
998
999 emin = mpfr_get_emin ();
1000 emax = mpfr_get_emax ();
1001 set_emin (MPFR_EMIN_MIN);
1002 set_emax (MPFR_EMAX_MAX);
1003
1004 mpq_init (qa);
1005 mpq_init (qb);
1006 mpfr_inits2 (32, x1, x2, y1, y2, y3, (mpfr_ptr) 0);
1007
1008 mpq_set_ui (qa, 3, 17);
1009 mpq_set_ui (qb, 17, 3);
1010 inex = mpfr_set_ui (x1, 7, MPFR_RNDN);
1011 MPFR_ASSERTN (inex == 0);
1012
1013 e = MPFR_EMAX_MAX - 3;
1014 inex = mpfr_set_ui_2exp (x2, 7, e, MPFR_RNDN); /* x2 = x1 * 2^e */
1015 MPFR_ASSERTN (inex == 0);
1016
1017 RND_LOOP(rnd)
1018 {
1019 mpfr_mul_q (y1, x1, qa, (mpfr_rnd_t) rnd);
1020 mpfr_div_q (y3, x1, qb, (mpfr_rnd_t) rnd);
1021 MPFR_ASSERTN (mpfr_equal_p (y1, y3));
1022 inex = mpfr_set_ui_2exp (y3, 1, e, MPFR_RNDN);
1023 MPFR_ASSERTN (inex == 0);
1024 inex = mpfr_mul (y3, y3, y1, MPFR_RNDN); /* y3 = y1 * 2^e */
1025 MPFR_ASSERTN (inex == 0);
1026 mpfr_mul_q (y2, x2, qa, (mpfr_rnd_t) rnd);
1027 if (! mpfr_equal_p (y2, y3))
1028 {
1029 printf ("Error 1 in bug_mul_div_q_20100818 (rnd = %d)\n", rnd);
1030 printf ("Expected "); mpfr_dump (y3);
1031 printf ("Got "); mpfr_dump (y2);
1032 exit (1);
1033 }
1034 mpfr_div_q (y2, x2, qb, (mpfr_rnd_t) rnd);
1035 if (! mpfr_equal_p (y2, y3))
1036 {
1037 printf ("Error 2 in bug_mul_div_q_20100818 (rnd = %d)\n", rnd);
1038 printf ("Expected "); mpfr_dump (y3);
1039 printf ("Got "); mpfr_dump (y2);
1040 exit (1);
1041 }
1042 }
1043
1044 e = MPFR_EMIN_MIN;
1045 inex = mpfr_set_ui_2exp (x2, 7, e, MPFR_RNDN); /* x2 = x1 * 2^e */
1046 MPFR_ASSERTN (inex == 0);
1047
1048 RND_LOOP(rnd)
1049 {
1050 mpfr_div_q (y1, x1, qa, (mpfr_rnd_t) rnd);
1051 mpfr_mul_q (y3, x1, qb, (mpfr_rnd_t) rnd);
1052 MPFR_ASSERTN (mpfr_equal_p (y1, y3));
1053 inex = mpfr_set_ui_2exp (y3, 1, e, MPFR_RNDN);
1054 MPFR_ASSERTN (inex == 0);
1055 inex = mpfr_mul (y3, y3, y1, MPFR_RNDN); /* y3 = y1 * 2^e */
1056 MPFR_ASSERTN (inex == 0);
1057 mpfr_div_q (y2, x2, qa, (mpfr_rnd_t) rnd);
1058 if (! mpfr_equal_p (y2, y3))
1059 {
1060 printf ("Error 3 in bug_mul_div_q_20100818 (rnd = %d)\n", rnd);
1061 printf ("Expected "); mpfr_dump (y3);
1062 printf ("Got "); mpfr_dump (y2);
1063 exit (1);
1064 }
1065 mpfr_mul_q (y2, x2, qb, (mpfr_rnd_t) rnd);
1066 if (! mpfr_equal_p (y2, y3))
1067 {
1068 printf ("Error 4 in bug_mul_div_q_20100818 (rnd = %d)\n", rnd);
1069 printf ("Expected "); mpfr_dump (y3);
1070 printf ("Got "); mpfr_dump (y2);
1071 exit (1);
1072 }
1073 }
1074
1075 mpq_clear (qa);
1076 mpq_clear (qb);
1077 mpfr_clears (x1, x2, y1, y2, y3, (mpfr_ptr) 0);
1078
1079 set_emin (emin);
1080 set_emax (emax);
1081 }
1082
1083 static void
reduced_expo_range(void)1084 reduced_expo_range (void)
1085 {
1086 mpfr_t x;
1087 mpz_t z;
1088 mpq_t q;
1089 mpfr_exp_t emin;
1090 int inex;
1091
1092 emin = mpfr_get_emin ();
1093 set_emin (4);
1094
1095 mpfr_init2 (x, 32);
1096
1097 mpz_init (z);
1098 mpfr_clear_flags ();
1099 inex = mpfr_set_ui (x, 17, MPFR_RNDN);
1100 MPFR_ASSERTN (inex == 0);
1101 mpz_set_ui (z, 3);
1102 inex = mpfr_mul_z (x, x, z, MPFR_RNDN);
1103 if (inex != 0 || MPFR_IS_NAN (x) || mpfr_cmp_ui (x, 51) != 0)
1104 {
1105 printf ("Error 1 in reduce_expo_range: expected 51 with inex = 0,"
1106 " got\n");
1107 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
1108 printf ("with inex = %d\n", inex);
1109 exit (1);
1110 }
1111 inex = mpfr_div_z (x, x, z, MPFR_RNDN);
1112 if (inex != 0 || MPFR_IS_NAN (x) || mpfr_cmp_ui (x, 17) != 0)
1113 {
1114 printf ("Error 2 in reduce_expo_range: expected 17 with inex = 0,"
1115 " got\n");
1116 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
1117 printf ("with inex = %d\n", inex);
1118 exit (1);
1119 }
1120 inex = mpfr_add_z (x, x, z, MPFR_RNDN);
1121 if (inex != 0 || MPFR_IS_NAN (x) || mpfr_cmp_ui (x, 20) != 0)
1122 {
1123 printf ("Error 3 in reduce_expo_range: expected 20 with inex = 0,"
1124 " got\n");
1125 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
1126 printf ("with inex = %d\n", inex);
1127 exit (1);
1128 }
1129 inex = mpfr_sub_z (x, x, z, MPFR_RNDN);
1130 if (inex != 0 || MPFR_IS_NAN (x) || mpfr_cmp_ui (x, 17) != 0)
1131 {
1132 printf ("Error 4 in reduce_expo_range: expected 17 with inex = 0,"
1133 " got\n");
1134 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
1135 printf ("with inex = %d\n", inex);
1136 exit (1);
1137 }
1138 MPFR_ASSERTN (__gmpfr_flags == 0);
1139 if (mpfr_cmp_z (x, z) <= 0)
1140 {
1141 printf ("Error 5 in reduce_expo_range: expected a positive value.\n");
1142 exit (1);
1143 }
1144 mpz_clear (z);
1145
1146 mpq_init (q);
1147 mpq_set_ui (q, 1, 1);
1148 mpfr_set_ui (x, 16, MPFR_RNDN);
1149 inex = mpfr_add_q (x, x, q, MPFR_RNDN);
1150 if (inex != 0 || MPFR_IS_NAN (x) || mpfr_cmp_ui (x, 17) != 0)
1151 {
1152 printf ("Error in reduce_expo_range for 16 + 1/1,"
1153 " got inex = %d and\nx = ", inex);
1154 mpfr_dump (x);
1155 exit (1);
1156 }
1157 inex = mpfr_sub_q (x, x, q, MPFR_RNDN);
1158 if (inex != 0 || MPFR_IS_NAN (x) || mpfr_cmp_ui (x, 16) != 0)
1159 {
1160 printf ("Error in reduce_expo_range for 17 - 1/1,"
1161 " got inex = %d and\nx = ", inex);
1162 mpfr_dump (x);
1163 exit (1);
1164 }
1165 mpq_clear (q);
1166
1167 mpfr_clear (x);
1168
1169 set_emin (emin);
1170 }
1171
1172 static void
addsubq_overflow_aux(mpfr_exp_t e)1173 addsubq_overflow_aux (mpfr_exp_t e)
1174 {
1175 mpfr_t x, y;
1176 mpq_t q;
1177 mpfr_exp_t emax;
1178 int inex;
1179 int rnd;
1180 int sign, sub;
1181
1182 MPFR_ASSERTN (e <= LONG_MAX);
1183 emax = mpfr_get_emax ();
1184 set_emax (e);
1185 mpfr_inits2 (16, x, y, (mpfr_ptr) 0);
1186 mpq_init (q);
1187
1188 mpfr_set_inf (x, 1);
1189 mpfr_nextbelow (x);
1190 mpq_set_ui (q, 1, 1);
1191
1192 for (sign = 0; sign <= 1; sign++)
1193 {
1194 for (sub = 0; sub <= 1; sub++)
1195 {
1196 RND_LOOP(rnd)
1197 {
1198 unsigned int flags, ex_flags;
1199 int inf;
1200
1201 inf = rnd == MPFR_RNDA ||
1202 rnd == (sign ? MPFR_RNDD : MPFR_RNDU);
1203 ex_flags = MPFR_FLAGS_INEXACT | (inf ? MPFR_FLAGS_OVERFLOW : 0);
1204 mpfr_clear_flags ();
1205 inex = sub ?
1206 mpfr_sub_q (y, x, q, (mpfr_rnd_t) rnd) :
1207 mpfr_add_q (y, x, q, (mpfr_rnd_t) rnd);
1208 flags = __gmpfr_flags;
1209 if (inex == 0 || flags != ex_flags ||
1210 (inf ? ! mpfr_inf_p (y) : ! mpfr_equal_p (x, y)))
1211 {
1212 printf ("Error in addsubq_overflow_aux(%ld),"
1213 " sign = %d, %s\n", (long) e, sign,
1214 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
1215 printf ("Got inex = %d, y = ", inex);
1216 mpfr_dump (y);
1217 printf ("Expected flags:");
1218 flags_out (ex_flags);
1219 printf ("Got flags: ");
1220 flags_out (flags);
1221 exit (1);
1222 }
1223 }
1224 mpq_neg (q, q);
1225 }
1226 mpfr_neg (x, x, MPFR_RNDN);
1227 mpq_neg (q, q);
1228 }
1229
1230 mpq_clear (q);
1231 mpfr_clears (x, y, (mpfr_ptr) 0);
1232 set_emax (emax);
1233 }
1234
1235 static void
addsubq_overflow(void)1236 addsubq_overflow (void)
1237 {
1238 addsubq_overflow_aux (4913);
1239 addsubq_overflow_aux (MPFR_EMAX_MAX);
1240 }
1241
1242 static void
coverage_mpfr_mul_q_20110218(void)1243 coverage_mpfr_mul_q_20110218 (void)
1244 {
1245 mpfr_t cmp, res, op1;
1246 mpq_t op2;
1247 int status;
1248
1249 mpfr_init2 (cmp, MPFR_PREC_MIN);
1250 mpfr_init2 (res, MPFR_PREC_MIN);
1251 mpfr_init_set_si (op1, 1, MPFR_RNDN);
1252
1253 mpq_init (op2);
1254 mpq_set_si (op2, 0, 0);
1255 mpz_set_si (mpq_denref (op2), 0);
1256
1257 status = mpfr_mul_q (res, op1, op2, MPFR_RNDN);
1258
1259 if ((status != 0) || (mpfr_cmp (cmp, res) != 0))
1260 {
1261 printf ("Results differ %d.\nres=", status);
1262 mpfr_dump (res);
1263 printf ("cmp=");
1264 mpfr_dump (cmp);
1265 exit (1);
1266 }
1267
1268 mpfr_set_si (op1, 1, MPFR_RNDN);
1269 mpq_set_si (op2, -1, 0);
1270
1271 status = mpfr_mul_q (res, op1, op2, MPFR_RNDN);
1272
1273 mpfr_set_inf (cmp, -1);
1274 if ((status != 0) || (mpfr_cmp(res, cmp) != 0))
1275 {
1276 printf ("mpfr_mul_q 1 * (-1/0) returned a wrong value:\n");
1277 printf (" expected ");
1278 mpfr_dump (cmp);
1279 printf (" got ");
1280 mpfr_dump (res);
1281 printf (" ternary value is %d\n", status);
1282 exit (1);
1283 }
1284
1285 mpq_clear (op2);
1286 mpfr_clear (op1);
1287 mpfr_clear (res);
1288 mpfr_clear (cmp);
1289 }
1290
1291 static void
coverage(void)1292 coverage (void)
1293 {
1294 mpfr_exp_t emax, emin;
1295 mpz_t z;
1296 mpfr_t x;
1297 int cmp;
1298
1299 mpz_init (z);
1300 mpfr_init2 (x, 5);
1301
1302 /* coverage for mpfr_cmp_z in case of overflow */
1303 emax = mpfr_get_emax ();
1304 set_emax (63);
1305 mpz_set_str (z, "9223372036854775808", 10); /* 2^63 */
1306 mpfr_set_ui_2exp (x, 1, mpfr_get_emax (), MPFR_RNDZ);
1307 /* x = (1-2^(-p))*2^emax */
1308 mpfr_clear_flags ();
1309 cmp = mpfr_cmp_z (x, z);
1310 MPFR_ASSERTN(cmp < 0);
1311 MPFR_ASSERTN(!mpfr_overflow_p ());
1312 set_emax (emax);
1313
1314 /* coverage for mpfr_cmp_z in case of underflow */
1315 mpz_set_str (z, "18446744073709551615", 10); /* 2^64-1 */
1316 emin = mpfr_get_emin ();
1317 set_emin (65); /* xmin = 2^64 */
1318 mpfr_set_ui_2exp (x, 1, 64, MPFR_RNDN);
1319 mpfr_clear_flags ();
1320 cmp = mpfr_cmp_z (x, z);
1321 MPFR_ASSERTN(cmp > 0);
1322 MPFR_ASSERTN(!mpfr_underflow_p ());
1323 set_emin (emin);
1324
1325 mpfr_clear (x);
1326 mpz_clear (z);
1327 }
1328
1329 int
main(int argc,char * argv[])1330 main (int argc, char *argv[])
1331 {
1332 tests_start_mpfr ();
1333
1334 coverage ();
1335 special ();
1336
1337 test_specialz (mpfr_add_z, mpz_add, "add");
1338 test_specialz (mpfr_sub_z, mpz_sub, "sub");
1339 test_specialz (mpfr_mul_z, mpz_mul, "mul");
1340 test_genericz (MPFR_PREC_MIN, 100, 100, mpfr_add_z, "add");
1341 test_genericz (MPFR_PREC_MIN, 100, 100, mpfr_sub_z, "sub");
1342 test_genericz (MPFR_PREC_MIN, 100, 100, mpfr_mul_z, "mul");
1343 test_genericz (MPFR_PREC_MIN, 100, 100, mpfr_div_z, "div");
1344 test_special2z (mpfr_z_sub, mpz_sub, "sub");
1345 test_generic2z (MPFR_PREC_MIN, 100, 100, mpfr_z_sub, "sub");
1346
1347 test_genericq (MPFR_PREC_MIN, 100, 100, mpfr_add_q, "add");
1348 test_genericq (MPFR_PREC_MIN, 100, 100, mpfr_sub_q, "sub");
1349 test_genericq (MPFR_PREC_MIN, 100, 100, mpfr_mul_q, "mul");
1350 test_genericq (MPFR_PREC_MIN, 100, 100, mpfr_div_q, "div");
1351 test_specialq (MPFR_PREC_MIN, 100, 100, mpfr_mul_q, mpq_mul, "mul");
1352 test_specialq (MPFR_PREC_MIN, 100, 100, mpfr_div_q, mpq_div, "div");
1353 test_specialq (MPFR_PREC_MIN, 100, 100, mpfr_add_q, mpq_add, "add");
1354 test_specialq (MPFR_PREC_MIN, 100, 100, mpfr_sub_q, mpq_sub, "sub");
1355
1356 test_cmp_z (MPFR_PREC_MIN, 100, 100);
1357 test_cmp_q (MPFR_PREC_MIN, 100, 100);
1358 test_cmp_f (MPFR_PREC_MIN, 100, 100);
1359
1360 check_for_zero ();
1361
1362 bug_mul_q_20100810 ();
1363 bug_div_q_20100810 ();
1364 bug_mul_div_q_20100818 ();
1365 reduced_expo_range ();
1366 addsubq_overflow ();
1367
1368 coverage_mpfr_mul_q_20110218 ();
1369
1370 tests_end_mpfr ();
1371 return 0;
1372 }
1373
1374 #else
1375
1376 int
main(void)1377 main (void)
1378 {
1379 return 77;
1380 }
1381
1382 #endif
1383