1284345Ssjg/* Test file for mpfr_urandom 2284345Ssjg 3284345SsjgCopyright 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4284345SsjgContributed by the Arenaire and Cacao projects, INRIA. 5284345Ssjg 6284345SsjgThis file is part of the GNU MPFR Library. 7284345Ssjg 8284345SsjgThe GNU MPFR Library is free software; you can redistribute it and/or modify 9284345Ssjgit under the terms of the GNU Lesser General Public License as published by 10284345Ssjgthe Free Software Foundation; either version 3 of the License, or (at your 11284345Ssjgoption) 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