1/* Test file for mpfr_urandomb 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_urandomb (long nbtests, mpfr_prec_t prec, int verbose) 27{ 28 mpfr_t x; 29 int *tab, size_tab, k, sh, xn; 30 double d, av = 0, var = 0, chi2 = 0, th; 31 mpfr_exp_t emin; 32 33 size_tab = (nbtests >= 1000 ? nbtests / 50 : 20); 34 tab = (int *) tests_allocate (size_tab * sizeof (int)); 35 for (k = 0; k < size_tab; k++) 36 tab[k] = 0; 37 38 mpfr_init2 (x, prec); 39 xn = 1 + (prec - 1) / mp_bits_per_limb; 40 sh = xn * mp_bits_per_limb - prec; 41 42 for (k = 0; k < nbtests; k++) 43 { 44 mpfr_urandomb (x, RANDS); 45 MPFR_ASSERTN (MPFR_IS_FP (x)); 46 /* check that lower bits are zero */ 47 if (MPFR_NOTZERO(x) && (MPFR_MANT(x)[0] & MPFR_LIMB_MASK(sh))) 48 { 49 printf ("Error: mpfr_urandomb() returns invalid numbers:\n"); 50 mpfr_dump (x); 51 exit (1); 52 } 53 d = mpfr_get_d1 (x); av += d; var += d*d; 54 tab[(int)(size_tab * d)]++; 55 } 56 57 /* coverage test */ 58 emin = mpfr_get_emin (); 59 set_emin (1); /* the generated number in [0,1[ is not in the exponent 60 range, except if it is zero */ 61 k = mpfr_urandomb (x, RANDS); 62 if (MPFR_IS_ZERO(x) == 0 && (k == 0 || mpfr_nan_p (x) == 0)) 63 { 64 printf ("Error in mpfr_urandomb, expected NaN, got "); 65 mpfr_dump (x); 66 exit (1); 67 } 68 set_emin (emin); 69 70 mpfr_clear (x); 71 72 if (verbose) 73 { 74 av /= nbtests; 75 var = (var / nbtests) - av * av; 76 77 th = (double) nbtests / size_tab; 78 printf ("Average = %.5f\nVariance = %.5f\n", av, var); 79 printf ("Repartition for urandomb. Each integer should be close to" 80 " %d.\n", (int) th); 81 82 for (k = 0; k < size_tab; k++) 83 { 84 chi2 += (tab[k] - th) * (tab[k] - th) / th; 85 printf ("%d ", tab[k]); 86 if (((unsigned int) (k+1) & 7) == 0) 87 printf ("\n"); 88 } 89 90 printf ("\nChi2 statistics value (with %d degrees of freedom) :" 91 " %.5f\n\n", size_tab - 1, chi2); 92 } 93 94 tests_free (tab, size_tab * sizeof (int)); 95 return; 96} 97 98#ifndef MPFR_USE_MINI_GMP 99/* Problem reported by Carl Witty: check mpfr_urandomb give similar results 100 on 32-bit and 64-bit machines. 101 We assume the default GMP random generator does not depend on the machine 102 word size, not on the GMP version. 103*/ 104static void 105bug20100914 (void) 106{ 107 mpfr_t x; 108 gmp_randstate_t s; 109 110#if __MPFR_GMP(4,2,0) 111# define C1 "0.895943" 112# define C2 "0.848824" 113#else 114# define C1 "0.479652" 115# define C2 "0.648529" 116#endif 117 118 gmp_randinit_default (s); 119 gmp_randseed_ui (s, 42); 120 mpfr_init2 (x, 17); 121 mpfr_urandomb (x, s); 122 if (mpfr_cmp_str1 (x, C1) != 0) 123 { 124 printf ("Error in bug20100914, expected " C1 ", got "); 125 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); 126 printf ("\n"); 127 exit (1); 128 } 129 mpfr_urandomb (x, s); 130 if (mpfr_cmp_str1 (x, C2) != 0) 131 { 132 printf ("Error in bug20100914, expected " C2 ", got "); 133 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); 134 printf ("\n"); 135 exit (1); 136 } 137 mpfr_clear (x); 138 gmp_randclear (s); 139} 140#endif 141 142int 143main (int argc, char *argv[]) 144{ 145 long nbtests; 146 mpfr_prec_t prec; 147 int verbose = 0; 148 149 tests_start_mpfr (); 150 151 if (argc > 1) 152 verbose = 1; 153 154 nbtests = 10000; 155 if (argc > 1) 156 { 157 long a = atol (argv[1]); 158 if (a != 0) 159 nbtests = a; 160 } 161 162 if (argc <= 2) 163 prec = 1000; 164 else 165 prec = atol (argv[2]); 166 167 test_urandomb (nbtests, prec, verbose); 168 169 if (argc == 1) /* check also small precision */ 170 { 171 test_urandomb (nbtests, 2, 0); 172 } 173 174#ifndef MPFR_USE_MINI_GMP 175 176 /* Since these tests assume a deterministic random generator, and this 177 is not implemented in mini-gmp, we omit them with mini-gmp. */ 178 179 bug20100914 (); 180 181#if __MPFR_GMP(4,2,0) 182 /* Get a non-zero fixed-point number whose first 32 bits are 0 with the 183 default GMP PRNG. This corresponds to the case cnt == 0 && k != 0 in 184 src/urandomb.c (fixed in r8762) with the 32-bit ABI. */ 185 { 186 gmp_randstate_t s; 187 mpfr_t x; 188 const char *str = "0.1010111100000000000000000000000000000000E-32"; 189 int k; 190 191 gmp_randinit_default (s); 192 gmp_randseed_ui (s, 4518); 193 mpfr_init2 (x, 40); 194 195 for (k = 0; k < 575123; k++) 196 { 197 mpfr_urandomb (x, s); 198 MPFR_ASSERTN (MPFR_IS_FP (x)); 199 } 200 201 if (mpfr_cmp_str (x, str, 2, MPFR_RNDN) != 0) 202 { 203 printf ("Error in test_urandomb:\n"); 204 printf ("Expected %s\n", str); 205 printf ("Got "); 206 mpfr_dump (x); 207 exit (1); 208 } 209 210 mpfr_clear (x); 211 gmp_randclear (s); 212 } 213#endif 214 215#endif 216 217 tests_end_mpfr (); 218 return 0; 219} 220