/* random_deviate routines for mpfr_erandom and mpfr_nrandom.
Copyright 2013-2018 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
http://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 mprf_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 runtime, 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;
}
}
}
/*
* return index [-1..127] of highest bit set. Return -1 if x = 0, 2 if x = 4,
* etc. (From Algorithms for programmers by Joerg Arndt.)
*/
static int
highest_bit_idx_alt (unsigned long x)
{
int r = 0;
if (x == 0)
return -1;
MPFR_ASSERTN (sizeof (unsigned long) * CHAR_BIT <= 128);
if (sizeof (unsigned long) * CHAR_BIT > 64)
{
/* handle 128-bit unsigned longs avoiding compiler warnings */
unsigned long y = x >> 16; y >>= 24; y >>= 24;
if (y) { 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;
}
/*
* return index [-1..63] of highest bit set.
* Return -1 if x = 0, 63 is if x = ~0 (for 64-bit unsigned long).
* See highest_bit_idx_alt too.
*/
static int
highest_bit_idx (unsigned long x)
{
/* this test should be evaluated at compile time */
if (sizeof (mp_limb_t) >= sizeof (unsigned long))
{
int cnt;
if (x == 0)
return -1;
count_leading_zeros (cnt, (mp_limb_t) x);
MPFR_ASSERTD (cnt <= GMP_NUMB_BITS - 1);
return GMP_NUMB_BITS - 1 - cnt;
}
else
return highest_bit_idx_alt (x);
}
/* 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;
}