1 /* Test file for mpfr_urandomb
2
3 Copyright 1999-2004, 2006-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba 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 https://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 "mpfr-test.h"
24
25 static void
26 test_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 */
104 static void
105 bug20100914 (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
142 int
143 main (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 }