1 /* Test file for mpfr_random_deviate
2
3 Copyright 2011-2023 Free Software Foundation, Inc.
4 Contributed by Charles Karney <charles@karney.com>, SRI International.
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 #include "random_deviate.h"
26
27 #define W 32 /* Must match value in random_deviate.c */
28
29 /* set random deviate rop from op */
30 static void
31 mpfr_random_deviate_set (mpfr_random_deviate_t rop, mpfr_random_deviate_t op)
32 {
33 rop->e = op->e;
34 rop->h = op->h;
35 mpz_set (rop->f, op->f);
36 }
37
38 /* set random deviate to fract * 2^-expt. expt must be a multiple
39 * of W and cannot be 0. fract must be in [0,2^W) */
40 static void
41 mpfr_random_deviate_ldexp (mpfr_random_deviate_t rop,
42 unsigned long fract, unsigned long expt)
43 {
44 rop->h = (expt > W ? 0ul : fract);
45 mpz_set_ui (rop->f, expt > W ? fract : 0ul);
46 rop->e = expt;
47 }
48
49 /* Test mpfr_random_deviate_less. With two initially equal deviates this
50 * should return true half the time. In order to execute additional code
51 * paths, the two deviates are repeatedly set equal and the test repeated (with
52 * now a longer fraction and with the test now triggering the sampling of an
53 * additional chunk. */
54 static void
55 test_compare (long nbtests, int verbose)
56 {
57 mpfr_random_deviate_t u, v;
58 int k, i, t1, t2;
59 long count;
60
61 mpfr_random_deviate_init (u);
62 mpfr_random_deviate_init (v);
63
64 count = 0;
65 for (k = 0; k < nbtests; ++k)
66 {
67 mpfr_random_deviate_reset (u);
68 mpfr_random_deviate_reset (v);
69 for (i = 0; i < 10; ++i)
70 {
71 t1 = mpfr_random_deviate_less (u, v, RANDS);
72 t2 = mpfr_random_deviate_less (u, v, RANDS);
73 if (t1 != t2)
74 {
75 printf ("Error: mpfr_random_deviate_less() inconsistent.\n");
76 exit (1);
77 }
78 if (t1)
79 ++count;
80 /* Force the test to sample an additional chunk */
81 mpfr_random_deviate_set (u, v);
82 }
83 t1 = mpfr_random_deviate_less (u, u, RANDS);
84 if (t1)
85 {
86 printf ("Error: mpfr_random_deviate_less() gives u < u.\n");
87 exit (1);
88 }
89 t1 = mpfr_random_deviate_tstbit (u, 0, RANDS);
90 if (t1)
91 {
92 printf ("Error: mpfr_random_deviate_tstbit() says 1 bit is on.\n");
93 exit (1);
94 }
95 }
96 mpfr_random_deviate_clear (v);
97 mpfr_random_deviate_clear (u);
98 if (verbose)
99 printf ("Fraction of true random_deviate_less = %.4f"
100 " (should be about 0.5)\n",
101 count / (double) (10 * nbtests));
102 }
103
104 /* A random_deviate should consistently return the same value at a given
105 * precision, even if intervening operations have caused the fraction to be
106 * extended. */
107 static void
108 test_value_consistency (long nbtests)
109 {
110 mpfr_t a1, a2, a3, b1, b2, b3;
111 mpfr_random_deviate_t u;
112 int inexa1, inexa2, inexa3, inexb1, inexb2, inexb3;
113 mpfr_prec_t prec1, prec2, prec3;
114 mpfr_rnd_t rnd;
115 long i;
116 unsigned n;
117 int neg;
118
119 /* Pick prec{1,2,3} random in [2,101] */
120 prec1 = 2 + gmp_urandomm_ui (RANDS, 100);
121 prec2 = 2 + gmp_urandomm_ui (RANDS, 100);
122 prec3 = 2 + gmp_urandomm_ui (RANDS, 100);
123
124 rnd = MPFR_RNDN;
125 mpfr_random_deviate_init (u);
126 mpfr_init2 (a1, prec1);
127 mpfr_init2 (b1, prec1);
128 mpfr_init2 (a2, prec2);
129 mpfr_init2 (b2, prec2);
130 mpfr_init2 (a3, prec3);
131 mpfr_init2 (b3, prec3);
132
133 for (i = 0; i < nbtests; ++i)
134 {
135 mpfr_random_deviate_reset (u);
136 n = gmp_urandomb_ui (RANDS, 4);
137 neg = gmp_urandomb_ui (RANDS, 1);
138 inexa1 = mpfr_random_deviate_value (neg, n, u, a1, RANDS, rnd);
139 inexa2 = mpfr_random_deviate_value (neg, n, u, a2, RANDS, rnd);
140 inexa3 = mpfr_random_deviate_value (neg, n, u, a3, RANDS, rnd);
141 inexb1 = mpfr_random_deviate_value (neg, n, u, b1, RANDS, rnd);
142 inexb2 = mpfr_random_deviate_value (neg, n, u, b2, RANDS, rnd);
143 inexb3 = mpfr_random_deviate_value (neg, n, u, b3, RANDS, rnd);
144 /* Of course a1, a2, and a3 should all be nearly equal. But more
145 * crucially (and easier to test), we need a1 == b1, etc. (This is not a
146 * trivial issue because the detailed mpfr operations giving b1 will be
147 * different than for a1, if, e.g., prec2 > prec1. */
148 if ( !( inexa1 == inexb1 && mpfr_equal_p (a1, b1) &&
149 inexa2 == inexb2 && mpfr_equal_p (a2, b2) &&
150 inexa3 == inexb3 && mpfr_equal_p (a3, b3) ) )
151 {
152 printf ("Error: random_deviate values are inconsistent.\n");
153 exit (1);
154 }
155 }
156 mpfr_random_deviate_clear (u);
157 mpfr_clears (a1, a2, a3, b1, b2, b3, (mpfr_ptr) 0);
158 }
159
160 /* Check that the values from random_deviate with different rounding modes are
161 * consistent. */
162 static void
163 test_value_round (long nbtests, mpfr_prec_t prec)
164 {
165 mpfr_t xn, xd, xu, xz, xa, t;
166 mpfr_random_deviate_t u;
167 int inexn, inexd, inexu, inexz, inexa, inext;
168 long i;
169 unsigned n;
170 int neg, s;
171
172 mpfr_random_deviate_init (u);
173 mpfr_inits2 (prec, xn, xd, xu, xz, xa, t, (mpfr_ptr) 0);
174
175 for (i = 0; i < nbtests; ++i)
176 {
177 mpfr_random_deviate_reset (u);
178 n = gmp_urandomb_ui (RANDS, 4);
179 neg = gmp_urandomb_ui (RANDS, 1);
180 s = neg ? -1 : 1;
181 inexn = mpfr_random_deviate_value (neg, n, u, xn, RANDS, MPFR_RNDN);
182 inexd = mpfr_random_deviate_value (neg, n, u, xd, RANDS, MPFR_RNDD);
183 inexu = mpfr_random_deviate_value (neg, n, u, xu, RANDS, MPFR_RNDU);
184 inexz = mpfr_random_deviate_value (neg, n, u, xz, RANDS, MPFR_RNDZ);
185 inexa = mpfr_random_deviate_value (neg, n, u, xa, RANDS, MPFR_RNDA);
186 inext = mpfr_set (t, xn, MPFR_RNDN);
187 /* Check inexact values */
188 if ( !( inexn != 0 && inext == 0 &&
189 inexd < 0 && inexu > 0 &&
190 inexz * s < 0 && inexa * s > 0 ) )
191 {
192 printf ("n %d t %d d %d u %d z %d a %d s %d\n",
193 inexn, inext, inexd, inexu, inexz, inexa, s);
194 printf ("Error: random_deviate has wrong values for inexact.\n");
195 exit (1);
196 }
197 if (inexn < 0)
198 mpfr_nextabove (t);
199 else
200 mpfr_nextbelow (t);
201 /* Check that x{d,u,z,a} == xn is the inexact flags match, else
202 * x{d,u,z,a} == t */
203 if ( !( mpfr_equal_p(xd, SAME_SIGN(inexn, inexd) ? xn : t) &&
204 mpfr_equal_p(xu, SAME_SIGN(inexn, inexu) ? xn : t) &&
205 mpfr_equal_p(xz, SAME_SIGN(inexn, inexz) ? xn : t) &&
206 mpfr_equal_p(xa, SAME_SIGN(inexn, inexa) ? xn : t) ) )
207 {
208 printf ("n %d t %d d %d u %d z %d a %d s %d\n",
209 inexn, inext, inexd, inexu, inexz, inexa, s);
210 printf ("n %.4f t %.4f d %.4f u %.4f z %.4f a %.4f\n",
211 mpfr_get_d (xn, MPFR_RNDN), mpfr_get_d (t, MPFR_RNDN),
212 mpfr_get_d (xd, MPFR_RNDN), mpfr_get_d (xu, MPFR_RNDN),
213 mpfr_get_d (xz, MPFR_RNDN), mpfr_get_d (xa, MPFR_RNDN));
214 printf ("Error: random_deviate rounds inconsistently.\n");
215 exit (1);
216 }
217 }
218 mpfr_random_deviate_clear (u);
219 mpfr_clears (xn, xd, xu, xz, xa, t, (mpfr_ptr) 0);
220 }
221
222 /* Test mpfr_random_deviate_value. Check for the leading bit in the number in
223 * various positions. */
224 static void
225 test_value (long nbtests, mpfr_prec_t prec, mpfr_rnd_t rnd,
226 int verbose)
227 {
228 mpfr_t x;
229 mpfr_random_deviate_t u;
230 int inexact, exact;
231 int i, k, b, neg;
232 unsigned long e, f, n;
233 long count, sum;
234
235 mpfr_random_deviate_init (u);
236 mpfr_init2 (x, prec);
237
238 count = 0; sum = 0;
239 exact = 0;
240
241 for (k = 0; k < nbtests; ++k)
242 {
243 for (i = 0; i < 32; ++i)
244 {
245 b = gmp_urandomm_ui (RANDS, 32) + 1; /* bits to sample in integer */
246 n = gmp_urandomb_ui (RANDS, b);
247 neg = gmp_urandomb_ui (RANDS, 1);
248 inexact = mpfr_random_deviate_value (neg, n, u, x, RANDS, rnd);
249 if (!inexact)
250 exact = 1;
251 if (inexact > 0)
252 ++count;
253 ++sum;
254 }
255 for (i = 0; i < 32; ++i)
256 {
257 b = gmp_urandomm_ui (RANDS, W) + 1; /* bits to sample in fraction */
258 f = gmp_urandomb_ui (RANDS, b);
259 e = W * (gmp_urandomm_ui (RANDS, 3) + 1);
260 mpfr_random_deviate_ldexp (u, f, e);
261 neg = gmp_urandomb_ui (RANDS, 1);
262 inexact = mpfr_random_deviate_value (neg, 0, u, x, RANDS, rnd);
263 if (!inexact)
264 exact = 1;
265 if (inexact > 0)
266 ++count;
267 ++sum;
268 }
269 if (exact)
270 {
271 printf ("Error: random_deviate() returns a zero ternary value.\n");
272 exit (1);
273 }
274 mpfr_random_deviate_reset (u);
275 }
276 mpfr_random_deviate_clear (u);
277 mpfr_clear (x);
278
279 if (verbose)
280 {
281 printf ("Fraction of inexact > 0 = %.4f", count / (double) (sum));
282 if (rnd == MPFR_RNDD)
283 printf (" should be exactly 0\n");
284 else if (rnd == MPFR_RNDU)
285 printf (" should be exactly 1\n");
286 else
287 printf (" should be about 0.5\n");
288 }
289 }
290
291 int
292 main (int argc, char *argv[])
293 {
294 long nbtests;
295 int verbose;
296 long a;
297
298 tests_start_mpfr ();
299
300 verbose = 0;
301 nbtests = 10;
302 if (argc > 1)
303 {
304 a = atol (argv[1]);
305 verbose = 1;
306 if (a != 0)
307 nbtests = a;
308 }
309
310 test_compare (nbtests, verbose);
311 test_value_consistency (nbtests);
312 test_value_round (nbtests, 2);
313 test_value_round (nbtests, 64);
314 test_value (nbtests, 2, MPFR_RNDD, verbose);
315 test_value (nbtests, 5, MPFR_RNDU, verbose);
316 test_value (nbtests, 24, MPFR_RNDN, verbose);
317 test_value (nbtests, 53, MPFR_RNDZ, verbose);
318 test_value (nbtests, 64, MPFR_RNDA, verbose);
319
320 tests_end_mpfr ();
321 return 0;
322 }