1/* Test file for mpfr_urandomb
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
28
29static void
30test_urandomb (long nbtests, mpfr_prec_t prec, 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
37  size_tab = (nbtests >= 1000 ? nbtests / 50 : 20);
38  tab = (int *) calloc (size_tab, sizeof(int));
39  if (tab == NULL)
40    {
41      fprintf (stderr, "trandom: can't allocate memory in test_urandomb\n");
42      exit (1);
43    }
44
45  mpfr_init2 (x, prec);
46  xn = 1 + (prec - 1) / mp_bits_per_limb;
47  sh = xn * mp_bits_per_limb - prec;
48
49  for (k = 0; k < nbtests; k++)
50    {
51      mpfr_urandomb (x, RANDS);
52      /* check that lower bits are zero */
53      if (MPFR_MANT(x)[0] & MPFR_LIMB_MASK(sh))
54        {
55          printf ("Error: mpfr_urandomb() returns invalid numbers:\n");
56          mpfr_print_binary (x); puts ("");
57          exit (1);
58        }
59      d = mpfr_get_d1 (x); av += d; var += d*d;
60      tab[(int)(size_tab * d)]++;
61    }
62
63  /* coverage test */
64  emin = mpfr_get_emin ();
65  set_emin (1); /* the generated number in [0,1[ is not in the exponent
66                        range, except if it is zero */
67  k = mpfr_urandomb (x, RANDS);
68  if (MPFR_IS_ZERO(x) == 0 && (k == 0 || mpfr_nan_p (x) == 0))
69    {
70      printf ("Error in mpfr_urandomb, expected NaN, got ");
71      mpfr_dump (x);
72      exit (1);
73    }
74  set_emin (emin);
75
76  mpfr_clear (x);
77  if (!verbose)
78    {
79      free(tab);
80      return;
81    }
82
83  av /= nbtests;
84  var = (var / nbtests) - av * av;
85
86  th = (double)nbtests / size_tab;
87  printf("Average = %.5f\nVariance = %.5f\n", av, var);
88  printf("Repartition for urandomb. Each integer should be close to %d.\n",
89         (int)th);
90
91  for (k = 0; k < size_tab; k++)
92    {
93      chi2 += (tab[k] - th) * (tab[k] - th) / th;
94      printf("%d ", tab[k]);
95      if (((k+1) & 7) == 0)
96        printf("\n");
97    }
98
99  printf("\nChi2 statistics value (with %d degrees of freedom) : %.5f\n\n",
100         size_tab - 1, chi2);
101
102  free(tab);
103  return;
104}
105
106int
107main (int argc, char *argv[])
108{
109  long nbtests;
110  mpfr_prec_t prec;
111  int verbose = 0;
112
113  tests_start_mpfr ();
114
115  if (argc > 1)
116    verbose = 1;
117
118  nbtests = 10000;
119  if (argc > 1)
120    {
121      long a = atol(argv[1]);
122      if (a != 0)
123        nbtests = a;
124    }
125
126  if (argc <= 2)
127    prec = 1000;
128  else
129    prec = atol(argv[2]);
130
131  test_urandomb (nbtests, prec, verbose);
132
133  if (argc == 1)  /* check also small precision */
134    {
135      test_urandomb (nbtests, 2, 0);
136    }
137
138  tests_end_mpfr ();
139  return 0;
140}
141