diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/random_deviate.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/random_deviate.c | 383 |
1 files changed, 0 insertions, 383 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/random_deviate.c b/Build/source/libs/mpfr/mpfr-src/src/random_deviate.c deleted file mode 100644 index 72531e774fa..00000000000 --- a/Build/source/libs/mpfr/mpfr-src/src/random_deviate.c +++ /dev/null @@ -1,383 +0,0 @@ -/* random_deviate routines for mpfr_erandom and mpfr_nrandom. - -Copyright 2013-2020 Free Software Foundation, Inc. -Contributed by Charles Karney <charles@karney.com>, 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_t x) -{ - mpz_init (x->f); - x->e = 0; -} - -/* reset to (0,1) */ -void -mpfr_random_deviate_reset (mpfr_random_deviate_t x) -{ - x->e = 0; -} - -/* deallocate */ -void -mpfr_random_deviate_clear (mpfr_random_deviate_t x) -{ - mpz_clear (x->f); -} - -/* swap two random deviates */ -void -mpfr_random_deviate_swap (mpfr_random_deviate_t x, mpfr_random_deviate_t 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_t 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_t 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_t 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_t x, mpfr_random_deviate_t 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_t x, mpfr_t 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; - - 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); - /* Is -x->e representable as a mpfr_exp_t? */ - MPFR_ASSERTN (x->e <= (mpfr_uexp_t)(-1) >> 1); - /* - * 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, -x->e, rnd); - mpz_clear (t); - return inex; -} |