1/* Test file for mpfr_urandomb
2
3Copyright 1999-2004, 2006-2023 Free Software Foundation, Inc.
4Contributed by the AriC and Caramba 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
20https://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 "mpfr-test.h"
24
25static void
26test_urandomb (long nbtests, mpfr_prec_t prec, int verbose)
27{
28  mpfr_t x;
29  int *tab, size_tab, k, sh, xn;
30  double d, av = 0, var = 0, chi2 = 0, th;
31  mpfr_exp_t emin;
32
33  size_tab = (nbtests >= 1000 ? nbtests / 50 : 20);
34  tab = (int *) tests_allocate (size_tab * sizeof (int));
35  for (k = 0; k < size_tab; k++)
36    tab[k] = 0;
37
38  mpfr_init2 (x, prec);
39  xn = 1 + (prec - 1) / mp_bits_per_limb;
40  sh = xn * mp_bits_per_limb - prec;
41
42  for (k = 0; k < nbtests; k++)
43    {
44      mpfr_urandomb (x, RANDS);
45      MPFR_ASSERTN (MPFR_IS_FP (x));
46      /* check that lower bits are zero */
47      if (MPFR_NOTZERO(x) && (MPFR_MANT(x)[0] & MPFR_LIMB_MASK(sh)))
48        {
49          printf ("Error: mpfr_urandomb() returns invalid numbers:\n");
50          mpfr_dump (x);
51          exit (1);
52        }
53      d = mpfr_get_d1 (x); av += d; var += d*d;
54      tab[(int)(size_tab * d)]++;
55    }
56
57  /* coverage test */
58  emin = mpfr_get_emin ();
59  set_emin (1); /* the generated number in [0,1[ is not in the exponent
60                        range, except if it is zero */
61  k = mpfr_urandomb (x, RANDS);
62  if (MPFR_IS_ZERO(x) == 0 && (k == 0 || mpfr_nan_p (x) == 0))
63    {
64      printf ("Error in mpfr_urandomb, expected NaN, got ");
65      mpfr_dump (x);
66      exit (1);
67    }
68  set_emin (emin);
69
70  mpfr_clear (x);
71
72  if (verbose)
73    {
74      av /= nbtests;
75      var = (var / nbtests) - av * av;
76
77      th = (double) nbtests / size_tab;
78      printf ("Average = %.5f\nVariance = %.5f\n", av, var);
79      printf ("Repartition for urandomb. Each integer should be close to"
80              " %d.\n", (int) th);
81
82      for (k = 0; k < size_tab; k++)
83        {
84          chi2 += (tab[k] - th) * (tab[k] - th) / th;
85          printf ("%d ", tab[k]);
86          if (((unsigned int) (k+1) & 7) == 0)
87            printf ("\n");
88        }
89
90      printf ("\nChi2 statistics value (with %d degrees of freedom) :"
91              " %.5f\n\n", size_tab - 1, chi2);
92    }
93
94  tests_free (tab, size_tab * sizeof (int));
95  return;
96}
97
98#ifndef MPFR_USE_MINI_GMP
99/* Problem reported by Carl Witty: check mpfr_urandomb give similar results
100   on 32-bit and 64-bit machines.
101   We assume the default GMP random generator does not depend on the machine
102   word size, not on the GMP version.
103*/
104static void
105bug20100914 (void)
106{
107  mpfr_t x;
108  gmp_randstate_t s;
109
110#if __MPFR_GMP(4,2,0)
111# define C1 "0.895943"
112# define C2 "0.848824"
113#else
114# define C1 "0.479652"
115# define C2 "0.648529"
116#endif
117
118  gmp_randinit_default (s);
119  gmp_randseed_ui (s, 42);
120  mpfr_init2 (x, 17);
121  mpfr_urandomb (x, s);
122  if (mpfr_cmp_str1 (x, C1) != 0)
123    {
124      printf ("Error in bug20100914, expected " C1 ", got ");
125      mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
126      printf ("\n");
127      exit (1);
128    }
129  mpfr_urandomb (x, s);
130  if (mpfr_cmp_str1 (x, C2) != 0)
131    {
132      printf ("Error in bug20100914, expected " C2 ", got ");
133      mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
134      printf ("\n");
135      exit (1);
136    }
137  mpfr_clear (x);
138  gmp_randclear (s);
139}
140#endif
141
142int
143main (int argc, char *argv[])
144{
145  long nbtests;
146  mpfr_prec_t prec;
147  int verbose = 0;
148
149  tests_start_mpfr ();
150
151  if (argc > 1)
152    verbose = 1;
153
154  nbtests = 10000;
155  if (argc > 1)
156    {
157      long a = atol (argv[1]);
158      if (a != 0)
159        nbtests = a;
160    }
161
162  if (argc <= 2)
163    prec = 1000;
164  else
165    prec = atol (argv[2]);
166
167  test_urandomb (nbtests, prec, verbose);
168
169  if (argc == 1)  /* check also small precision */
170    {
171      test_urandomb (nbtests, 2, 0);
172    }
173
174#ifndef MPFR_USE_MINI_GMP
175
176  /* Since these tests assume a deterministic random generator, and this
177     is not implemented in mini-gmp, we omit them with mini-gmp. */
178
179  bug20100914 ();
180
181#if __MPFR_GMP(4,2,0)
182  /* Get a non-zero fixed-point number whose first 32 bits are 0 with the
183     default GMP PRNG. This corresponds to the case cnt == 0 && k != 0 in
184     src/urandomb.c (fixed in r8762) with the 32-bit ABI. */
185  {
186    gmp_randstate_t s;
187    mpfr_t x;
188    const char *str = "0.1010111100000000000000000000000000000000E-32";
189    int k;
190
191    gmp_randinit_default (s);
192    gmp_randseed_ui (s, 4518);
193    mpfr_init2 (x, 40);
194
195    for (k = 0; k < 575123; k++)
196      {
197        mpfr_urandomb (x, s);
198        MPFR_ASSERTN (MPFR_IS_FP (x));
199      }
200
201    if (mpfr_cmp_str (x, str, 2, MPFR_RNDN) != 0)
202      {
203        printf ("Error in test_urandomb:\n");
204        printf ("Expected %s\n", str);
205        printf ("Got      ");
206        mpfr_dump (x);
207        exit (1);
208      }
209
210    mpfr_clear (x);
211    gmp_randclear (s);
212  }
213#endif
214
215#endif
216
217  tests_end_mpfr ();
218  return 0;
219}
220