1/* Test file for mpfr_urandom 2 3Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao 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 20http://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 <stdio.h> 24#include <stdlib.h> 25 26#include "mpfr-test.h" 27 28static void 29test_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 157int 158main (int argc, char *argv[]) 159{ 160 long nbtests; 161 mpfr_prec_t prec; 162 int verbose = 0; 163 int rnd; 164 long bit_index; 165 166 tests_start_mpfr (); 167 168 if (argc > 1) 169 verbose = 1; 170 171 nbtests = 10000; 172 if (argc > 1) 173 { 174 long a = atol(argv[1]); 175 if (a != 0) 176 nbtests = a; 177 } 178 179 if (argc <= 2) 180 prec = 1000; 181 else 182 prec = atol(argv[2]); 183 184 if (argc <= 3) 185 bit_index = -1; 186 else 187 { 188 bit_index = atol(argv[3]); 189 if (bit_index >= prec) 190 { 191 printf ("Warning. Cannot compute the bit frequency: the given bit " 192 "index (= %ld) is not less than the precision (= %ld).\n", 193 bit_index, prec); 194 bit_index = -1; 195 } 196 } 197 198 RND_LOOP(rnd) 199 { 200 test_urandom (nbtests, prec, (mpfr_rnd_t) rnd, bit_index, verbose); 201 202 if (argc == 1) /* check also small precision */ 203 { 204 test_urandom (nbtests, 2, (mpfr_rnd_t) rnd, -1, 0); 205 } 206 } 207 208 tests_end_mpfr (); 209 return 0; 210} 211