/* random_deviate routines for mpfr_erandom and mpfr_nrandom. Copyright 2013-2021 Free Software Foundation, Inc. Contributed by Charles Karney , SRI International. This file is part of the GNU MPFR Library. The GNU MPFR Library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. The GNU MPFR Library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ /* * A mpfr_random_deviate represents the initial portion e bits of a random * deviate uniformly distributed in (0,1) as * * typedef struct { * unsigned long e; // bits in the fraction * unsigned long h; // the high W bits of the fraction * mpz_t f; // the rest of the fraction * } mpfr_random_deviate_t[1]; * * e is always a multiple of RANDOM_CHUNK. The first RANDOM_CHUNK bits, the * high fraction, are held in an unsigned long, h, and the rest are held in an * mpz_t, f. The data in h is undefined if e == 0 and, similarly the data in f * is undefined if e <= RANDOM_CHUNK. */ #define MPFR_NEED_LONGLONG_H #include "random_deviate.h" /* * RANDOM_CHUNK can be picked in the range 1 <= RANDOM_CHUNK <= 64. Low values * of RANDOM_CHUNK are good for testing, since they are more likely to make * bugs obvious. For portability, pick RANDOM_CHUNK <= 32 (since an unsigned * long may only hold 32 bits). For reproducibility across platforms, * standardize on RANDOM_CHUNK = 32. * * When RANDOM_CHUNK = 32, this representation largely avoids manipulating * mpz's (until the final cast to an mpfr is done). In addition * mpfr_random_deviate_less usually entails just a single comparison of * unsigned longs. In this way, we can stick with the published interface for * extracting portions of an mpz (namely through mpz_tstbit) without hurting * efficiency. */ #if !defined(RANDOM_CHUNK) /* note: for MPFR, we could use RANDOM_CHUNK = 32 or 64 according to the number of bits per limb, but we use 32 everywhere to get reproducible results on 32-bit and 64-bit computers */ #define RANDOM_CHUNK 32 /* Require 1 <= RANDOM_CHUNK <= 32; recommend 32 */ #endif #define W RANDOM_CHUNK /* W is just an shorter name for RANDOM_CHUNK */ /* allocate and set to (0,1) */ void mpfr_random_deviate_init (mpfr_random_deviate_ptr x) { mpz_init (x->f); x->e = 0; } /* reset to (0,1) */ void mpfr_random_deviate_reset (mpfr_random_deviate_ptr x) { x->e = 0; } /* deallocate */ void mpfr_random_deviate_clear (mpfr_random_deviate_ptr x) { mpz_clear (x->f); } /* swap two random deviates */ void mpfr_random_deviate_swap (mpfr_random_deviate_ptr x, mpfr_random_deviate_ptr y) { mpfr_random_size_t s; unsigned long t; /* swap x->e and y->e */ s = x->e; x->e = y->e; y->e = s; /* swap x->h and y->h */ t = x->h; x->h = y->h; y->h = t; /* swap x->f and y->f */ mpz_swap (x->f, y->f); } /* ensure x has at least k bits */ static void random_deviate_generate (mpfr_random_deviate_ptr x, mpfr_random_size_t k, gmp_randstate_t r, mpz_t t) { /* Various compile time checks on mpfr_random_deviate_t */ /* Check that the h field of a mpfr_random_deviate_t can hold W bits */ MPFR_STAT_STATIC_ASSERT (W > 0 && W <= sizeof (unsigned long) * CHAR_BIT); /* Check mpfr_random_size_t can hold 32 bits and a mpfr_uprec_t. This * ensures that max(mpfr_random_size_t) exceeds MPFR_PREC_MAX by at least * 2^31 because mpfr_prec_t is a signed version of mpfr_uprec_t. This allows * random deviates with many leading zeros in the fraction to be handled * correctly. */ MPFR_STAT_STATIC_ASSERT (sizeof (mpfr_random_size_t) * CHAR_BIT >= 32 && sizeof (mpfr_random_size_t) >= sizeof (mpfr_uprec_t)); /* Finally, at run time, check that k is not too big. e is set to ceil(k/W)*W * and we require that this allows x->e + 1 in random_deviate_leading_bit to * be computed without overflow. */ MPFR_ASSERTN (k <= (mpfr_random_size_t)(-((int) W + 1))); /* if t is non-null, it is used as a temporary */ if (x->e >= k) return; if (x->e == 0) { x->h = gmp_urandomb_ui (r, W); /* Generate the high fraction */ x->e = W; if (x->e >= k) return; /* Maybe that's it? */ } if (t) { /* passed a mpz_t so compute needed bits in one call to mpz_urandomb */ k = ((k + (W-1)) / W) * W; /* Round up to multiple of W */ k -= x->e; /* The number of new bits */ mpz_urandomb (x->e == W ? x->f : t, r, k); /* Copy directly to x->f? */ if (x->e > W) { mpz_mul_2exp (x->f, x->f, k); mpz_add (x->f, x->f, t); } x->e += k; } else { /* no mpz_t so compute the bits W at a time via gmp_urandomb_ui */ while (x->e < k) { unsigned long w = gmp_urandomb_ui (r, W); if (x->e == W) mpz_set_ui (x->f, w); else { mpz_mul_2exp (x->f, x->f, W); mpz_add_ui (x->f, x->f, w); } x->e += W; } } } #ifndef MPFR_LONG_WITHIN_LIMB /* a long does not fit in a mp_limb_t */ /* * return index [0..127] of highest bit set. Return 0 if x = 1, 2 if x = 4, * etc. Assume x > 0. (From Algorithms for programmers by Joerg Arndt.) */ static int highest_bit_idx (unsigned long x) { unsigned long y; int r = 0; MPFR_ASSERTD(x > 0); MPFR_STAT_STATIC_ASSERT (sizeof (unsigned long) * CHAR_BIT <= 128); /* A compiler with VRP (like GCC) will optimize and not generate any code for the following lines if unsigned long has at most 64 values bits. */ y = ((x >> 16) >> 24) >> 24; /* portable x >> 64 */ if (y != 0) { x = y; r += 64; } if (x & ~0xffffffffUL) { x >>= 16; x >>= 16; r +=32; } if (x & 0xffff0000UL) { x >>= 16; r += 16; } if (x & 0x0000ff00UL) { x >>= 8; r += 8; } if (x & 0x000000f0UL) { x >>= 4; r += 4; } if (x & 0x0000000cUL) { x >>= 2; r += 2; } if (x & 0x00000002UL) { r += 1; } return r; } #else /* a long fits in a mp_limb_t */ /* * return index [0..63] of highest bit set. Assume x > 0. * Return 0 if x = 1, 63 is if x = ~0 (for 64-bit unsigned long). * See alternate code above too. */ static int highest_bit_idx (unsigned long x) { int cnt; MPFR_ASSERTD(x > 0); count_leading_zeros (cnt, (mp_limb_t) x); MPFR_ASSERTD (cnt <= GMP_NUMB_BITS - 1); return GMP_NUMB_BITS - 1 - cnt; } #endif /* MPFR_LONG_WITHIN_LIMB */ /* return position of leading bit, counting from 1 */ static mpfr_random_size_t random_deviate_leading_bit (mpfr_random_deviate_ptr x, gmp_randstate_t r) { mpfr_random_size_t l; random_deviate_generate (x, W, r, 0); if (x->h) return W - highest_bit_idx (x->h); random_deviate_generate (x, 2 * W, r, 0); while (mpz_sgn (x->f) == 0) random_deviate_generate (x, x->e + 1, r, 0); l = x->e + 1 - mpz_sizeinbase (x->f, 2); /* Guard against a ridiculously long string of leading zeros in the fraction; * probability of this happening is 2^(-2^31). In particular ensure that * p + 1 + l in mpfr_random_deviate_value doesn't overflow with p = * MPFR_PREC_MAX. */ MPFR_ASSERTN (l + 1 < (mpfr_random_size_t)(-MPFR_PREC_MAX)); return l; } /* return kth bit of fraction, representing 2^-k */ int mpfr_random_deviate_tstbit (mpfr_random_deviate_ptr x, mpfr_random_size_t k, gmp_randstate_t r) { if (k == 0) return 0; random_deviate_generate (x, k, r, 0); if (k <= W) return (x->h >> (W - k)) & 1UL; return mpz_tstbit (x->f, x->e - k); } /* compare two random deviates, x < y */ int mpfr_random_deviate_less (mpfr_random_deviate_ptr x, mpfr_random_deviate_ptr y, gmp_randstate_t r) { mpfr_random_size_t k = 1; if (x == y) return 0; random_deviate_generate (x, W, r, 0); random_deviate_generate (y, W, r, 0); if (x->h != y->h) return x->h < y->h; /* Compare the high fractions */ k += W; for (; ; ++k) { /* Compare the rest of the fraction bit by bit */ int a = mpfr_random_deviate_tstbit (x, k, r); int b = mpfr_random_deviate_tstbit (y, k, r); if (a != b) return a < b; } } /* set mpfr_t z = (neg ? -1 : 1) * (n + x) */ int mpfr_random_deviate_value (int neg, unsigned long n, mpfr_random_deviate_ptr x, mpfr_ptr z, gmp_randstate_t r, mpfr_rnd_t rnd) { /* r is used to add as many bits as necessary to match the precision of z */ int s; mpfr_random_size_t l; /* The leading bit is 2^(s*l) */ mpfr_random_size_t p = mpfr_get_prec (z); /* Number of bits in result */ mpz_t t; int inex; mpfr_exp_t negxe; if (n == 0) { s = -1; l = random_deviate_leading_bit (x, r); /* l > 0 */ } else { s = 1; l = highest_bit_idx (n); /* l >= 0 */ } /* * Leading bit is 2^(s*l); thus the trailing bit in result is 2^(s*l-p+1) = * 2^-(p-1-s*l). For the sake of illustration, take l = 0 and p = 4, thus * bits through the 1/8 position need to be generated; assume that these bits * are 1.010 = 10/8 which represents a deviate in the range (10,11)/8. * * If the rounding mode is one of RNDZ, RNDU, RNDD, RNDA, we add a 1 bit to * the result to give 1.0101 = (10+1/2)/8. When this is converted to a MPFR * the result is rounded to 10/8, 11/8, 10/8, 11/8, respectively, and the * inexact flag is set to -1, 1, -1, 1. * * If the rounding mode is RNDN, an additional random bit must be generated * to determine if the result is in (10,10+1/2)/8 or (10+1/2,11)/8. Assume * that this random bit is 0, so the result is 1.0100 = (10+0/2)/8. Then an * additional 1 bit is added to give 1.010101 = (10+1/4)/8. This last bit * avoids the "round ties to even rule" (because there are no ties) and sets * the inexact flag so that the result is 10/8 with the inexact flag = 1. * * Here we always generate at least 2 additional random bits, so that bit * position 2^-(p+1-s*l) is generated. (The result often contains more * random bits than this because random bits are added in batches of W and * because additional bits may have been required in the process of * generating the random deviate.) The integer and all the bits in the * fraction are then copied into an mpz, the least significant bit is * unconditionally set to 1, the sign is set, and the result together with * the exponent -x->e is used to generate an mpfr using mpfr_set_z_2exp. * * If random bits were very expensive, we would only need to generate to the * 2^-(p-1-s*l) bit (no extra bits) for the RNDZ, RNDU, RNDD, RNDA modes and * to the 2^-(p-s*l) bit (1 extra bit) for RNDN. By always generating 2 bits * we save on some bit shuffling when formed the mpz to be converted to an * mpfr. The implementation of the RandomNumber class in RandomLib * illustrates the more parsimonious approach (which was taken to allow * accurate counts of the number of random digits to be made). */ mpz_init (t); /* * This is the only call to random_deviate_generate where a mpz_t is passed * (because an arbitrarily large number of bits may need to be generated). */ if ((s > 0 && p + 1 > l) || (s < 0 && p + 1 + l > 0)) random_deviate_generate (x, s > 0 ? p + 1 - l : p + 1 + l, r, t); if (n == 0) { /* Since the minimum prec is 2 we know that x->h has been generated. */ mpz_set_ui (t, x->h); /* Set high fraction */ } else { mpz_set_ui (t, n); /* The integer part */ if (x->e > 0) { mpz_mul_2exp (t, t, W); /* Shift to allow for high fraction */ mpz_add_ui (t, t, x->h); /* Add high fraction */ } } if (x->e > W) { mpz_mul_2exp (t, t, x->e - W); /* Shift to allow for low fraction */ mpz_add (t, t, x->f); /* Add low fraction */ } /* * We could trim off any excess bits here by shifting rightward. This is an * unnecessary complication. */ mpz_setbit (t, 0); /* Set the trailing bit so result is always inexact */ if (neg) mpz_neg (t, t); /* Portable version of the negation of x->e, with a check of overflow. */ if (MPFR_UNLIKELY (x->e > MPFR_EXP_MAX)) { /* Overflow, except when x->e = MPFR_EXP_MAX + 1 = - MPFR_EXP_MIN. */ MPFR_ASSERTN (MPFR_EXP_MIN + MPFR_EXP_MAX == -1 && x->e == (mpfr_random_size_t) MPFR_EXP_MAX + 1); negxe = MPFR_EXP_MIN; } else negxe = - (mpfr_exp_t) x->e; /* * Let mpfr_set_z_2exp do all the work of rounding to the requested * precision, setting overflow/underflow flags, and returning the right * inexact value. */ inex = mpfr_set_z_2exp (z, t, negxe, rnd); mpz_clear (t); return inex; }