1/* Test file for mpfr_urandom 2 3Copyright 1999-2004, 2006-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba 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 20https://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 "mpfr-test.h" 24 25static void 26test_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 184static void 185underflow_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 279static void 280test_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 312static void 313overflow_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. */ 381static void 382bug20100914 (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 */ 420static void 421bug20170123 (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. */ 442static void 443reprod_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*/ 522static void 523reprod_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 602int 603main (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} 670