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