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 } 323