1/* Test file for mpfr_random_deviate
2
3Copyright 2011-2023 Free Software Foundation, Inc.
4Contributed by Charles Karney <charles@karney.com>, SRI International.
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
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 */
30static void
31mpfr_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) */
40static void
41mpfr_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. */
54static void
55test_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. */
107static void
108test_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. */
162static void
163test_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. */
224static void
225test_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
291int
292main (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}
323