tests/turandom.c

   1 /* Test file for mpfr_urandom
   2
   3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011, 2012 Free Software Foundation, Inc.
   4 Contributed by the AriC and Caramel 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 }