1 /* Test file for mpfr_urandom 2 3 Copyright 1999-2004, 2006-2016 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 http://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 <stdio.h> 24 #include <stdlib.h> 25 26 #include "mpfr-test.h" 27 28 static void 29 test_urandom (long nbtests, mpfr_prec_t prec, mpfr_rnd_t rnd, long bit_index, 30 int verbose) 31 { 32 mpfr_t x; 33 int *tab, size_tab, k, sh, xn; 34 double d, av = 0, var = 0, chi2 = 0, th; 35 mpfr_exp_t emin; 36 mp_size_t limb_index = 0; 37 mp_limb_t limb_mask = 0; 38 long count = 0; 39 int i; 40 int inex = 1; 41 42 size_tab = (nbtests >= 1000 ? nbtests / 50 : 20); 43 tab = (int *) calloc (size_tab, sizeof(int)); 44 if (tab == NULL) 45 { 46 fprintf (stderr, "trandom: can't allocate memory in test_urandom\n"); 47 exit (1); 48 } 49 50 mpfr_init2 (x, prec); 51 xn = 1 + (prec - 1) / mp_bits_per_limb; 52 sh = xn * mp_bits_per_limb - prec; 53 if (bit_index >= 0 && bit_index < prec) 54 { 55 /* compute the limb index and limb mask to fetch the bit #bit_index */ 56 limb_index = (prec - bit_index) / mp_bits_per_limb; 57 i = 1 + bit_index - (bit_index / mp_bits_per_limb) * mp_bits_per_limb; 58 limb_mask = MPFR_LIMB_ONE << (mp_bits_per_limb - i); 59 } 60 61 for (k = 0; k < nbtests; k++) 62 { 63 i = mpfr_urandom (x, RANDS, rnd); 64 inex = (i != 0) && inex; 65 /* check that lower bits are zero */ 66 if (MPFR_MANT(x)[0] & MPFR_LIMB_MASK(sh) && !MPFR_IS_ZERO (x)) 67 { 68 printf ("Error: mpfr_urandom() returns invalid numbers:\n"); 69 mpfr_print_binary (x); puts (""); 70 exit (1); 71 } 72 /* check that the value is in [0,1] */ 73 if (mpfr_cmp_ui (x, 0) < 0 || mpfr_cmp_ui (x, 1) > 0) 74 { 75 printf ("Error: mpfr_urandom() returns number outside [0, 1]:\n"); 76 mpfr_print_binary (x); puts (""); 77 exit (1); 78 } 79 80 d = mpfr_get_d1 (x); av += d; var += d*d; 81 i = (int)(size_tab * d); 82 if (d == 1.0) i --; 83 tab[i]++; 84 85 if (limb_mask && (MPFR_MANT (x)[limb_index] & limb_mask)) 86 count ++; 87 } 88 89 if (inex == 0) 90 { 91 /* one call in the loop pretended to return an exact number! */ 92 printf ("Error: mpfr_urandom() returns a zero ternary value.\n"); 93 exit (1); 94 } 95 96 /* coverage test */ 97 emin = mpfr_get_emin (); 98 for (k = 0; k < 5; k++) 99 { 100 set_emin (k+1); 101 inex = mpfr_urandom (x, RANDS, rnd); 102 if (( (rnd == MPFR_RNDZ || rnd == MPFR_RNDD) 103 && (!MPFR_IS_ZERO (x) || inex != -1)) 104 || ((rnd == MPFR_RNDU || rnd == MPFR_RNDA) 105 && (mpfr_cmp_ui (x, 1 << k) != 0 || inex != +1)) 106 || (rnd == MPFR_RNDN 107 && (k > 0 || mpfr_cmp_ui (x, 1 << k) != 0 || inex != +1) 108 && (!MPFR_IS_ZERO (x) || inex != -1))) 109 { 110 printf ("Error: mpfr_urandom() do not handle correctly a restricted" 111 " exponent range.\nrounding mode: %s\nternary value: %d\n" 112 "random value: ", mpfr_print_rnd_mode (rnd), inex); 113 mpfr_print_binary (x); puts (""); 114 exit (1); 115 } 116 } 117 set_emin (emin); 118 119 mpfr_clear (x); 120 if (!verbose) 121 { 122 free(tab); 123 return; 124 } 125 126 av /= nbtests; 127 var = (var / nbtests) - av * av; 128 129 th = (double)nbtests / size_tab; 130 printf ("Average = %.5f\nVariance = %.5f\n", av, var); 131 printf ("Repartition for urandom with rounding mode %s. " 132 "Each integer should be close to %d.\n", 133 mpfr_print_rnd_mode (rnd), (int)th); 134 135 for (k = 0; k < size_tab; k++) 136 { 137 chi2 += (tab[k] - th) * (tab[k] - th) / th; 138 printf("%d ", tab[k]); 139 if (((k+1) & 7) == 0) 140 printf("\n"); 141 } 142 143 printf("\nChi2 statistics value (with %d degrees of freedom) : %.5f\n", 144 size_tab - 1, chi2); 145 146 if (limb_mask) 147 printf ("Bit #%ld is set %ld/%ld = %.1f %% of time\n", 148 bit_index, count, nbtests, count * 100.0 / nbtests); 149 150 puts (""); 151 152 free(tab); 153 return; 154 } 155 156 /* problem reported by Carl Witty */ 157 static void 158 bug20100914 (void) 159 { 160 mpfr_t x; 161 gmp_randstate_t s; 162 163 #if __MPFR_GMP(4,2,0) 164 # define C1 "0.8488312" 165 # define C2 "0.8156509" 166 #else 167 # define C1 "0.6485367" 168 # define C2 "0.9362717" 169 #endif 170 171 gmp_randinit_default (s); 172 gmp_randseed_ui (s, 42); 173 mpfr_init2 (x, 17); 174 mpfr_urandom (x, s, MPFR_RNDN); 175 if (mpfr_cmp_str1 (x, C1) != 0) 176 { 177 printf ("Error in bug20100914, expected " C1 ", got "); 178 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); 179 printf ("\n"); 180 exit (1); 181 } 182 mpfr_urandom (x, s, MPFR_RNDN); 183 if (mpfr_cmp_str1 (x, C2) != 0) 184 { 185 printf ("Error in bug20100914, expected " C2 ", got "); 186 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); 187 printf ("\n"); 188 exit (1); 189 } 190 mpfr_clear (x); 191 gmp_randclear (s); 192 } 193 194 int 195 main (int argc, char *argv[]) 196 { 197 long nbtests; 198 mpfr_prec_t prec; 199 int verbose = 0; 200 int rnd; 201 long bit_index; 202 203 tests_start_mpfr (); 204 205 if (argc > 1) 206 verbose = 1; 207 208 nbtests = 10000; 209 if (argc > 1) 210 { 211 long a = atol(argv[1]); 212 if (a != 0) 213 nbtests = a; 214 } 215 216 if (argc <= 2) 217 prec = 1000; 218 else 219 prec = atol(argv[2]); 220 221 if (argc <= 3) 222 bit_index = -1; 223 else 224 { 225 bit_index = atol(argv[3]); 226 if (bit_index >= prec) 227 { 228 printf ("Warning. Cannot compute the bit frequency: the given bit " 229 "index (= %ld) is not less than the precision (= %ld).\n", 230 bit_index, (long) prec); 231 bit_index = -1; 232 } 233 } 234 235 RND_LOOP(rnd) 236 { 237 test_urandom (nbtests, prec, (mpfr_rnd_t) rnd, bit_index, verbose); 238 239 if (argc == 1) /* check also small precision */ 240 { 241 test_urandom (nbtests, 2, (mpfr_rnd_t) rnd, -1, 0); 242 } 243 } 244 245 bug20100914 (); 246 247 tests_end_mpfr (); 248 return 0; 249 } 250