summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/random_deviate.c
diff options
context:
space:
mode:
authorKarl Berry <karl@freefriends.org>2021-02-25 19:22:25 +0000
committerKarl Berry <karl@freefriends.org>2021-02-25 19:22:25 +0000
commitad547a6b5986815fda458221149728d9d9ab1d87 (patch)
tree16296910eb3eca724371474ea9aea3994dc69614 /Build/source/libs/mpfr/mpfr-src/src/random_deviate.c
parent947b43de3dd21d58ccc2ffadefc4441ea1c2a813 (diff)
restore Build,TODO from r57911
git-svn-id: svn://tug.org/texlive/trunk@57915 c570f23f-e606-0410-a88d-b1316a301751
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.c383
1 files changed, 383 insertions, 0 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
new file mode 100644
index 00000000000..72531e774fa
--- /dev/null
+++ b/Build/source/libs/mpfr/mpfr-src/src/random_deviate.c
@@ -0,0 +1,383 @@
+/* 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;
+}