1 /* Test file for mpfr_urandom
2
3 Copyright 1999-2004, 2006-2023 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 static void
26 test_urandom (long nbtests, mpfr_prec_t prec, mpfr_rnd_t rnd, long bit_index,
27 int verbose)
28 {
29 mpfr_t x;
30 int *tab, size_tab, k, sh, xn;
31 double d, av = 0, var = 0, chi2 = 0, th;
32 mpfr_exp_t emin;
33 mp_size_t limb_index = 0;
34 mp_limb_t limb_mask = 0;
35 long count = 0;
36 int i;
37 int inex = 1;
38 mpfr_flags_t ex_flags, flags;
39
40 size_tab = (nbtests >= 1000 ? nbtests / 50 : 20);
41 tab = (int *) tests_allocate (size_tab * sizeof (int));
42 for (k = 0; k < size_tab; k++)
43 tab[k] = 0;
44
45 mpfr_init2 (x, prec);
46 xn = 1 + (prec - 1) / mp_bits_per_limb;
47 sh = xn * mp_bits_per_limb - prec;
48 if (bit_index >= 0 && bit_index < prec)
49 {
50 /* compute the limb index and limb mask to fetch the bit #bit_index */
51 limb_index = (prec - bit_index) / mp_bits_per_limb;
52 i = 1 + bit_index - (bit_index / mp_bits_per_limb) * mp_bits_per_limb;
53 limb_mask = MPFR_LIMB_ONE << (mp_bits_per_limb - i);
54 }
55
56 for (k = 0; k < nbtests; k++)
57 {
58 mpfr_clear_flags ();
59 ex_flags = MPFR_FLAGS_INEXACT;
60 i = mpfr_urandom (x, RANDS, rnd);
61 flags = __gmpfr_flags;
62 inex = (i != 0) && inex;
63 /* check that lower bits are zero */
64 if (MPFR_MANT(x)[0] & MPFR_LIMB_MASK(sh) && !MPFR_IS_ZERO (x))
65 {
66 printf ("Error: mpfr_urandom() returns invalid numbers:\n");
67 mpfr_dump (x);
68 exit (1);
69 }
70 /* check that the value is in [0,1] */
71 if (mpfr_cmp_ui (x, 0) < 0 || mpfr_cmp_ui (x, 1) > 0)
72 {
73 printf ("Error: mpfr_urandom() returns number outside [0, 1]:\n");
74 mpfr_dump (x);
75 exit (1);
76 }
77 /* check the flags (an underflow is theoretically possible, but
78 impossible in practice due to the huge exponent range) */
79 if (flags != ex_flags)
80 {
81 printf ("Error: mpfr_urandom() returns incorrect flags.\n");
82 printf ("Expected ");
83 flags_out (ex_flags);
84 printf ("Got ");
85 flags_out (flags);
86 exit (1);
87 }
88
89 d = mpfr_get_d1 (x);
90 av += d;
91 var += d*d;
92 i = (int) (size_tab * d);
93 if (d == 1.0)
94 i--;
95 MPFR_ASSERTN (i < size_tab);
96 tab[i]++;
97
98 if (limb_mask && (MPFR_MANT (x)[limb_index] & limb_mask))
99 count ++;
100 }
101
102 if (inex == 0)
103 {
104 /* one call in the loop pretended to return an exact number! */
105 printf ("Error: mpfr_urandom() returns a zero ternary value.\n");
106 exit (1);
107 }
108
109 /* coverage test */
110 emin = mpfr_get_emin ();
111 for (k = 0; k < 5; k++)
112 {
113 set_emin (k+1);
114 ex_flags = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
115 for (i = 0; i < 5; i++)
116 {
117 mpfr_clear_flags ();
118 inex = mpfr_urandom (x, RANDS, rnd);
119 flags = __gmpfr_flags;
120 if (k > 0 && flags != ex_flags)
121 {
122 printf ("Error: mpfr_urandom() returns incorrect flags"
123 " for emin = %d (i = %d).\n", k+1, i);
124 printf ("Expected ");
125 flags_out (ex_flags);
126 printf ("Got ");
127 flags_out (flags);
128 exit (1);
129 }
130 if (( (rnd == MPFR_RNDZ || rnd == MPFR_RNDD)
131 && (!MPFR_IS_ZERO (x) || inex != -1))
132 || ((rnd == MPFR_RNDU || rnd == MPFR_RNDA)
133 && (mpfr_cmp_ui (x, 1 << k) != 0 || inex != +1))
134 || (rnd == MPFR_RNDN
135 && (k > 0 || mpfr_cmp_ui (x, 1 << k) != 0 || inex != +1)
136 && (!MPFR_IS_ZERO (x) || inex != -1)))
137 {
138 printf ("Error: mpfr_urandom() does not handle correctly"
139 " a restricted exponent range.\nemin = %d\n"
140 "rounding mode: %s\nternary value: %d\nrandom value: ",
141 k+1, mpfr_print_rnd_mode (rnd), inex);
142 mpfr_dump (x);
143 exit (1);
144 }
145 }
146 }
147 set_emin (emin);
148
149 mpfr_clear (x);
150
151 if (verbose)
152 {
153 av /= nbtests;
154 var = (var / nbtests) - av * av;
155
156 th = (double)nbtests / size_tab;
157 printf ("Average = %.5f\nVariance = %.5f\n", av, var);
158 printf ("Repartition for urandom with rounding mode %s. "
159 "Each integer should be close to %d.\n",
160 mpfr_print_rnd_mode (rnd), (int) th);
161
162 for (k = 0; k < size_tab; k++)
163 {
164 chi2 += (tab[k] - th) * (tab[k] - th) / th;
165 printf("%d ", tab[k]);
166 if (((unsigned int) (k+1) & 7) == 0)
167 printf("\n");
168 }
169
170 printf("\nChi2 statistics value (with %d degrees of freedom) : %.5f\n",
171 size_tab - 1, chi2);
172
173 if (limb_mask)
174 printf ("Bit #%ld is set %ld/%ld = %.1f %% of time\n",
175 bit_index, count, nbtests, count * 100.0 / nbtests);
176
177 puts ("");
178 }
179
180 tests_free (tab, size_tab * sizeof (int));
181 return;
182 }
183
184 static void
185 underflow_tests (void)
186 {
187 mpfr_t x;
188 mpfr_exp_t emin;
189 int i, k;
190 int rnd;
191
192 emin = mpfr_get_emin ();
193 mpfr_init2 (x, 4);
194
195 for (i = 2; i >= -4; i--)
196 RND_LOOP (rnd)
197 for (k = 0; k < 100; k++)
198 {
199 mpfr_flags_t ex_flags, flags;
200 int inex;
201
202 if (i >= 2)
203 {
204 /* Always underflow when emin >= 2, i.e. when the minimum
205 representable positive number is >= 2. */
206 ex_flags = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
207 }
208 else
209 {
210 #ifndef MPFR_USE_MINI_GMP
211 gmp_randstate_t s;
212
213 /* Since the unrounded random number does not depend on
214 the current exponent range, we can detect underflow
215 in a range larger than the one that will be tested. */
216 gmp_randinit_set (s, mpfr_rands);
217 mpfr_clear_flags ();
218 mpfr_urandom (x, s, (mpfr_rnd_t) rnd);
219 gmp_randclear (s);
220 ex_flags = MPFR_FLAGS_INEXACT;
221 if (MPFR_IS_ZERO (x) || mpfr_get_exp (x) < i)
222 ex_flags |= MPFR_FLAGS_UNDERFLOW;
223 #else
224 /* Do not test the flags. */
225 ex_flags = 0;
226 #endif
227 }
228
229 set_emin (i);
230 mpfr_clear_flags ();
231 inex = mpfr_urandom (x, mpfr_rands, (mpfr_rnd_t) rnd);
232 flags = __gmpfr_flags;
233 MPFR_ASSERTN (mpfr_inexflag_p ());
234 set_emin (emin);
235
236 if (MPFR_IS_NEG (x))
237 {
238 printf ("Error in underflow_tests: got a negative sign"
239 " for i=%d rnd=%s k=%d.\n",
240 i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k);
241 exit (1);
242 }
243
244 if (MPFR_IS_ZERO (x))
245 {
246 if (rnd == MPFR_RNDU || rnd == MPFR_RNDA)
247 {
248 printf ("Error in underflow_tests: the value cannot"
249 " be 0 for i=%d rnd=%s k=%d.\n",
250 i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k);
251 exit (1);
252 }
253 }
254
255 if (inex == 0 || (MPFR_IS_ZERO (x) && inex > 0))
256 {
257 printf ("Error in underflow_tests: incorrect inex (%d)"
258 " for i=%d rnd=%s k=%d.\n", inex,
259 i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k);
260 exit (1);
261 }
262
263 if (ex_flags != 0 && flags != ex_flags)
264 {
265 printf ("Error in underflow_tests: incorrect flags"
266 " for i=%d rnd=%s k=%d.\n",
267 i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k);
268 printf ("Expected ");
269 flags_out (ex_flags);
270 printf ("Got ");
271 flags_out (flags);
272 exit (1);
273 }
274 }
275
276 mpfr_clear (x);
277 }
278
279 static void
280 test_underflow (int verbose)
281 {
282 mpfr_t x;
283 mpfr_exp_t emin = mpfr_get_emin ();
284 long i, exp[6] = {0, 0, 0, 0, 0, 0};
285
286 mpfr_init2 (x, 2);
287 set_emin (-3);
288 #define N 1000000
289 for (i = 0; i < N; i++)
290 {
291 mpfr_urandom (x, RANDS, MPFR_RNDN);
292 if (mpfr_zero_p (x))
293 exp[5] ++;
294 else
295 /* exp=1 is possible if the generated number is 0.111111... */
296 exp[1-mpfr_get_exp(x)] ++;
297 }
298 if (verbose)
299 printf ("exp=1:%.3f(%.3f) 0:%.3f(%.3f) -1:%.3f(%.3f) -2:%.3f(%.3f) "
300 "-3:%.3f(%.3f) zero:%.3f(%.3f)\n",
301 100.0 * (double) exp[0] / (double) N, 12.5,
302 100.0 * (double) exp[1] / (double) N, 43.75,
303 100.0 * (double) exp[2] / (double) N, 21.875,
304 100.0 * (double) exp[3] / (double) N, 10.9375,
305 100.0 * (double) exp[4] / (double) N, 7.8125,
306 100.0 * (double) exp[5] / (double) N, 3.125);
307 mpfr_clear (x);
308 set_emin (emin);
309 #undef N
310 }
311
312 static void
313 overflow_tests (void)
314 {
315 mpfr_t x;
316 mpfr_exp_t emax;
317 int i, k;
318 int inex;
319 int rnd;
320 mpfr_flags_t ex_flags, flags;
321
322 emax = mpfr_get_emax ();
323 mpfr_init2 (x, 4);
324 ex_flags = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT; /* if overflow */
325 for (i = -4; i <= 0; i++)
326 {
327 set_emax (i);
328 RND_LOOP (rnd)
329 for (k = 0; k < 100; k++)
330 {
331 mpfr_clear_flags ();
332 inex = mpfr_urandom (x, mpfr_rands, (mpfr_rnd_t) rnd);
333 flags = __gmpfr_flags;
334 MPFR_ASSERTN (mpfr_inexflag_p ());
335 if (MPFR_IS_NEG (x))
336 {
337 printf ("Error in overflow_tests: got a negative sign"
338 " for i=%d rnd=%s k=%d.\n",
339 i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k);
340 exit (1);
341 }
342 if (MPFR_IS_INF (x))
343 {
344 if (rnd == MPFR_RNDD || rnd == MPFR_RNDZ)
345 {
346 printf ("Error in overflow_tests: the value cannot"
347 " be +inf for i=%d rnd=%s k=%d.\n",
348 i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k);
349 exit (1);
350 }
351 if (flags != ex_flags)
352 {
353 printf ("Error in overflow_tests: incorrect flags"
354 " for i=%d rnd=%s k=%d.\n",
355 i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k);
356 printf ("Expected ");
357 flags_out (ex_flags);
358 printf ("Got ");
359 flags_out (flags);
360 exit (1);
361 }
362 }
363 if (inex == 0 || (MPFR_IS_INF (x) && inex < 0))
364 {
365 printf ("Error in overflow_tests: incorrect inex (%d)"
366 " for i=%d rnd=%s k=%d.\n", inex,
367 i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k);
368 exit (1);
369 }
370 }
371 }
372 mpfr_clear (x);
373 set_emax (emax);
374 }
375
376 #ifndef MPFR_USE_MINI_GMP
377
378 /* Problem reported by Carl Witty. This test assumes the random generator
379 used by GMP is deterministic (for a given seed). We need to distinguish
380 two cases since the random generator changed in GMP 4.2.0. */
381 static void
382 bug20100914 (void)
383 {
384 mpfr_t x;
385 gmp_randstate_t s;
386
387 #if __MPFR_GMP(4,2,0)
388 # define C1 "0.8488312"
389 # define C2 "0.8156509"
390 #else
391 # define C1 "0.6485367"
392 # define C2 "0.9362717"
393 #endif
394
395 gmp_randinit_default (s);
396 gmp_randseed_ui (s, 42);
397 mpfr_init2 (x, 17);
398 mpfr_urandom (x, s, MPFR_RNDN);
399 if (mpfr_cmp_str1 (x, C1) != 0)
400 {
401 printf ("Error in bug20100914, expected " C1 ", got ");
402 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
403 printf ("\n");
404 exit (1);
405 }
406 mpfr_urandom (x, s, MPFR_RNDN);
407 if (mpfr_cmp_str1 (x, C2) != 0)
408 {
409 printf ("Error in bug20100914, expected " C2 ", got ");
410 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
411 printf ("\n");
412 exit (1);
413 }
414 mpfr_clear (x);
415 gmp_randclear (s);
416 }
417
418 /* non-regression test for bug reported by Trevor Spiteri
419 https://sympa.inria.fr/sympa/arc/mpfr/2017-01/msg00020.html */
420 static void
421 bug20170123 (void)
422 {
423 #if __MPFR_GMP(4,2,0)
424 mpfr_t x;
425 mpfr_exp_t emin;
426 gmp_randstate_t s;
427
428 emin = mpfr_get_emin ();
429 set_emin (-7);
430 mpfr_init2 (x, 53);
431 gmp_randinit_default (s);
432 gmp_randseed_ui (s, 398);
433 mpfr_urandom (x, s, MPFR_RNDN);
434 MPFR_ASSERTN(mpfr_cmp_ui_2exp (x, 1, -8) == 0);
435 mpfr_clear (x);
436 gmp_randclear (s);
437 set_emin (emin);
438 #endif
439 }
440
441 /* Reproducibility test with several rounding modes and exponent ranges. */
442 static void
443 reprod_rnd_exp (void)
444 {
445 int i;
446
447 for (i = 0; i < 10; i++)
448 {
449 gmp_randstate_t s1;
450 mpfr_prec_t prec;
451 mpfr_t x1, x2, y;
452 mp_limb_t v;
453 int r;
454
455 prec = MPFR_PREC_MIN + (randlimb () % 200);
456 mpfr_inits2 (prec, x1, x2, y, (mpfr_ptr) 0);
457
458 gmp_randinit_set (s1, mpfr_rands);
459 mpfr_urandom (x1, mpfr_rands, MPFR_RNDZ);
460 mpfr_rand_raw (&v, mpfr_rands, GMP_NUMB_BITS);
461 mpfr_set (x2, x1, MPFR_RNDN);
462 mpfr_nextabove (x2);
463 /* The real number is between x1 and x2. */
464
465 RND_LOOP (r)
466 {
467 gmp_randstate_t s2;
468 mpfr_rnd_t rr = (mpfr_rnd_t) r;
469 mp_limb_t w;
470 mpfr_ptr t[2];
471 int k, nk = 0;
472
473 gmp_randinit_set (s2, s1);
474 mpfr_urandom (y, s2, rr);
475 mpfr_rand_raw (&w, s2, GMP_NUMB_BITS);
476 if (w != v)
477 {
478 printf ("Error in reprod_rnd_exp for i=%d rnd=%s: different "
479 "PRNG state\n", i, mpfr_print_rnd_mode (rr));
480 exit (1);
481 }
482
483 if (! MPFR_IS_LIKE_RNDA (rr, 0))
484 t[nk++] = x1;
485 if (! MPFR_IS_LIKE_RNDZ (rr, 0))
486 t[nk++] = x2;
487 MPFR_ASSERTN (nk == 1 || nk == 2);
488
489 if (!(mpfr_equal_p (y, t[0]) || (nk > 1 && mpfr_equal_p (y, t[1]))))
490 {
491 printf ("Error in reprod_rnd_exp for i=%d rnd=%s:\n",
492 i, mpfr_print_rnd_mode (rr));
493 printf ("Expected%s\n", nk > 1 ? " one of" : "");
494 for (k = 0; k < nk; k++)
495 {
496 printf (" ");
497 mpfr_dump (t[k]);
498 }
499 printf ("Got\n ");
500 mpfr_dump (y);
501 exit (1);
502 }
503
504 gmp_randclear (s2);
505 }
506
507 mpfr_clears (x1, x2, y, (mpfr_ptr) 0);
508 gmp_randclear (s1);
509 }
510 }
511
512 /* Reproducibility test: check that the behavior does not depend on
513 the platform ABI or MPFR version (new, incompatible MPFR versions
514 may introduce changes, in which case the hardcoded values should
515 depend on MPFR_VERSION).
516 It is not necessary to test with different rounding modes and
517 exponent ranges as this has already been done in reprod_rnd_exp.
518 We do not need to check the status of the PRNG after mpfr_urandom
519 since this is done implicitly by comparing the next value, except
520 for the last itaration.
521 */
522 static void
523 reprod_abi (void)
524 {
525 #if __MPFR_GMP(4,2,0)
526 #define N 6
527 /* Run this program with the MPFR_REPROD_ABI_OUTPUT environment variable
528 set to get the array of strings. */
529 const char *t[5 * N] = {
530 "1.0@-1",
531 "3.0@-1",
532 "7.0@-1",
533 "9.0@-1",
534 "c.0@-1",
535 "4.385434c0@-1",
536 "1.9a018734@-1",
537 "8.26547780@-1",
538 "a.fd334198@-1",
539 "9.aa11d5f00@-1",
540 "d.aa9a32fd0a801ac0@-1",
541 "c.eb47074368ec6340@-1",
542 "9.7dbe2ced88ae4c30@-1",
543 "d.03218ea6704a42c0@-1",
544 "7.1530156aac800d980@-1",
545 "e.270121b1d74aea9029ccc740@-1",
546 "5.614fc2d9ca3917107609e2e0@-1",
547 "5.15417c51af272232181d6a40@-1",
548 "f.dfb431dd6533c004b6d3c590@-1",
549 "4.345f96fd2929d41eb278a4f40@-1",
550 "a.804590c6449cd8c83bae31f31f7a4100@-1",
551 "a.0a2b318d3c99911a45e4cf33847d3680@-1",
552 "2.89f6127c19092d7a1808b1842b296400@-1",
553 "2.1db4fc00348ca1531983fe4bd4cdf6d2@-1",
554 "5.2d90f11ed710425ebe549a95decbb6540@-1",
555 "8.ca35d1302cf369e03c2a58bf2ce5cff8307f0bc0@-1",
556 "3.a22bae632c32f2a7a67a1fa78a93f5e84f9caa40@-1",
557 "f.370a36febed972dbb47f2503f7e08a651edbf120@-1",
558 "d.0764d7a38c206eeba6ffe8cf39d777194f5c9200@-1",
559 "a.1a312f0bb16db20c4783c6438725ed5d6dff6af80@-1"
560 };
561 gmp_randstate_t s;
562 int generate, i;
563
564 /* We must hardcode the seed to be able to compare with hardcoded values. */
565 gmp_randinit_default (s);
566 gmp_randseed_ui (s, 17);
567
568 generate = getenv ("MPFR_REPROD_ABI_OUTPUT") != NULL;
569
570 for (i = 0; i < 5 * N; i++)
571 {
572 mpfr_prec_t prec;
573 mpfr_t x;
574
575 prec = i < 5 ? MPFR_PREC_MIN + i : (i / 5) * 32 + (i % 5) - 2;
576 mpfr_init2 (x, prec);
577 mpfr_urandom (x, s, MPFR_RNDN);
578 if (generate)
579 {
580 printf (" \"");
581 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDZ);
582 printf (i < 5 * N - 1 ? "\",\n" : "\"\n");
583 }
584 else if (mpfr_cmp_str (x, t[i], 16, MPFR_RNDN) != 0)
585 {
586 printf ("Error in reprod_abi for i=%d\n", i);
587 printf ("Expected %s\n", t[i]);
588 printf ("Got ");
589 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDZ);
590 printf ("\n");
591 exit (1);
592 }
593 mpfr_clear (x);
594 }
595
596 gmp_randclear (s);
597 #endif
598 }
599
600 #endif
601
602 int
603 main (int argc, char *argv[])
604 {
605 long nbtests;
606 mpfr_prec_t prec;
607 int verbose = 0;
608 int rnd;
609 long bit_index;
610
611 tests_start_mpfr ();
612
613 if (argc > 1)
614 verbose = 1;
615
616 nbtests = 10000;
617 if (argc > 1)
618 {
619 long a = atol(argv[1]);
620 if (a != 0)
621 nbtests = a;
622 }
623
624 if (argc <= 2)
625 prec = 1000;
626 else
627 prec = atol(argv[2]);
628
629 if (argc <= 3)
630 bit_index = -1;
631 else
632 {
633 bit_index = atol(argv[3]);
634 if (bit_index >= prec)
635 {
636 printf ("Warning. Cannot compute the bit frequency: the given bit "
637 "index (= %ld) is not less than the precision (= %ld).\n",
638 bit_index, (long) prec);
639 bit_index = -1;
640 }
641 }
642
643 RND_LOOP(rnd)
644 {
645 test_urandom (nbtests, prec, (mpfr_rnd_t) rnd, bit_index, verbose);
646
647 if (argc == 1) /* check also small precision */
648 {
649 test_urandom (nbtests, MPFR_PREC_MIN, (mpfr_rnd_t) rnd, -1, 0);
650 }
651 }
652
653 underflow_tests ();
654 overflow_tests ();
655
656 #ifndef MPFR_USE_MINI_GMP
657 /* Since these tests assume a deterministic random generator, and
658 this is not implemented in mini-gmp, we omit it with mini-gmp. */
659 bug20100914 ();
660 bug20170123 ();
661 reprod_rnd_exp ();
662 reprod_abi ();
663 #endif
664
665 test_underflow (verbose);
666
667 tests_end_mpfr ();
668 return 0;
669 }