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