/*-
* Copyright (c) 2018 Taylor R. Campbell
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
/*
* If you take the uniform distribution on [0,1] or (0,1), i.e. the
* Lebesgue measure, and round it to the nearest floating-point number,
* the distribution on the result is the natural notion of the uniform
* distribution on the floating-point numbers in [0,1].
*
* Note that even if you round a sample from the open interval (0,1) of
* real numbers, the result might be rounded to 0 or 1, if it is in (0,
* 2^{emin - p - 1}] or [1 - eps/4, 1):
*
* - The interval near 0 is of negligible size in any format larger
* than binary16 -- less than 2^-126 in binary32 even if you use
* flush-to-zero, less than 2^-1022 in binary64 -- so you don't need
* to worry about it and can safely compute, e.g., 1/x or log(x).
*
* - However, the probability of getting a real number in the interval
* near 1 is 2^-54 in binary64, which is small but not negligible, so
* I recommend not rejecting it unless you must avoid 1, for instance
* because you need to compute log1p(-x).
*
* In infinite-precision floating-point, the exponent is geometrically
* distributed in {-1, -2, -3, ...}, and the fractional part of the
* significand is a uniformly distributed infinite bit string. If we
* round this to finite precision, the exponent remains the same
* (unless rounded to subnormal or zero) and the significand is rounded
* to p bits. The outcome of rounding is not actually affected by more
* than a finite number of significand bits, with one tricky exception,
* so you don't need to sample an infinite stream of bits and round it;
* to get this, you can:
*
* 1. flip a coin until you get one heads, to pick a geometrically
* distributed exponent;
*
* 2. flip a coin >p+1 times to choose the uniformly distributed
* fractional part of a significand, and round it to p bits with a
* trick to avoid a bias toward `even'.
*
* If you drew p bits for the fractional part of the significand, there
* would be no rounding, and you would exclude from the possible
* significands the interval (2 - eps/2, 2] of real numbers that are
* rounded to 2 -- small but not nonnegligible. If you draw >p bits
* and round it, then no real number in [1, 2] is excluded.
*
* However, there is a slight bias toward `even' outcomes -- that is,
* outcomes whose least significant bit is zero -- because there's a
* nonzero probability that the finite significand you choose is a tie,
* exactly between two floating-point numbers, which is broken by
* choosing the even option. This event has probability zero in
* Lebesgue measure. You can simply break the ties by drawing >p+1
* bits and forcing the least significant bit to be 1.
*
* The attached file uniform01.c implements this algorithm in the
* function sample_uniform_01, with parameters chosen for IEEE 754
* binary64 (`double-precision') floating-point.
*
* As it happens, the variant algorithm that chooses exactly p bits and
* _doesn't_ round them gives a uniform distribution on [0, 1 - eps/4),
* which incidentally effectively never returns zero, and so yields
* what you might think of `the uniform distribution on floating-point
* numbers in (0, 1)'. I named it sample_uniform_01_open, and added
* some paranoia to clamp the result at the smallest subnormal in case
* your PRNG is broken. But it also works to just do rejection
* sampling on sample_uniform_01 and reduce the amount of code to worry
* about.
*
* Note: This code does not run in constant time -- the time it takes
* is proportional to the magnitude of the exponent. (In principle, it
* may also take variable time for subnormal arithmetic, but the
* probability of that is negligible for any format larger than
* binary16.)
*
* You could make it run in constant time by always flipping the
* maximum number of coins for an exponent, and counting the leading
* zero bits. To improve performance, you can safely truncate the
* exponent at, say, 128 -- if you flip 128 coins and you get tails
* every time, either your coin flipper is broken or there's a glitch
* in the matrix.
*/
#include
#include
#include
/**
* Draw an integer in [0, 2^32) uniformly at random.
*/
#include /* arc4random */
uint32_t
crypto_rand_uint32(void)
{
return arc4random(); /* XXX */
}
/**
* Count number of one bits in 32-bit word
*/
unsigned
bitcount32(uint32_t x)
{
/* Count two-bit groups. */
x -= (x >> 1) & UINT32_C(0x55555555);
/* Count four-bit groups. */
x = ((x >> 2) & UINT32_C(0x33333333)) + (x & UINT32_C(0x33333333));
/* Count eight-bit groups. */
x = (x + (x >> 4)) & UINT32_C(0x0f0f0f0f);
/* Sum all eight-bit groups, and extract the sum. */
return (x * UINT32_C(0x01010101)) >> 24;
}
/**
* Count leading zeros in 32-bit word
*/
unsigned
clz32(uint32_t x)
{
/* Round up to a power of two. */
x |= x >> 1;
x |= x >> 2;
x |= x >> 4;
x |= x >> 8;
x |= x >> 16;
/* Subtract count of one bits from 32. */
return (32 - bitcount32(x));
}
/*
* The following sample_uniform_01 is tailored for IEEE 754 binary64
* floating-point or smaller. It can be adapted to larger
* floating-point formats like i387 80-bit or IEEE 754 binary128, but
* it may require sampling more bits.
*/
CTASSERT(FLT_RADIX == 2);
CTASSERT(-DBL_MIN_EXP <= 1021);
CTASSERT(DBL_MANT_DIG <= 53);
/**
* Sample a floating-point number in [0, 1] with uniform distribution.
*
* Note that the probability of returning 0 is less than 2^-1074, so
* callers need not check for it. However, callers that cannot handle
* rounding to 1 must deal with that, because it occurs with
* probability 2^-54, which is small but nonnegligible.
*/
double
sample_uniform_01(void)
{
uint32_t z, x, hi, lo;
double s;
/*
* Draw an exponent, geometrically distributed, but give up if
* we get a run of more than 1088 zeros, which really means the
* system is broken.
*/
z = 0;
while ((x = crypto_rand_uint32()) == 0) {
if (z >= 1088)
/* Your bit sampler is broken. Go home. */
return 0;
z += 32;
}
z += clz32(x);
/*
* Pick 32-bit halves of an odd normalized significand.
* Picking it odd breaks ties in the subsequent rounding, which
* occur only with measure zero in the uniform distribution on
* [0, 1].
*/
hi = crypto_rand_uint32() | UINT32_C(0x80000000);
lo = crypto_rand_uint32() | UINT32_C(0x00000001);
/* Round to nearest scaled significand in [2^63, 2^64]. */
s = hi*(double)4294967296 + lo;
/* Rescale into [1/2, 1] and apply exponent in one swell foop. */
return s * ldexp(1, -(64 + z));
}
/*
* The following sample_uniform_01_open is tailored for IEEE 754
* binary64 floating-point exactly. It can be adapted to smaller or
* larger floating-point formats, but it must fit the bits exactly to
* avoid rounding.
*/
CTASSERT(FLT_RADIX == 2);
CTASSERT(-DBL_MIN_EXP == 1021);
CTASSERT(DBL_MANT_DIG == 53);
/**
* Sample a floating-point number in (2^{emin-p-1}, 1 - eps/4) with
* uniform distribution. The outcome lies in [2^{emin-p}, 1 - eps/2].
* This is sometimes described as `(0, 1)', but you are guaranteed not
* to get the floating-point rounding of any real number in the
* intervals [0, 2^{emin-p-1}) or (1 - eps/4, 1].
*
* Use this only if you are sure you can't handle the value 1, because
* the width of the interval that is rounded to 1, 2^-54 in binary64,
* is small but not negligible. sample_uniform_01() is already
* effectively guaranteed not to return 0.
*/
double
sample_uniform_01_open(void)
{
#if 0
/* Rejection sampling. */
double x;
/* 0 is effectively impossible -- probability <2^-1074. */
do {
x = uniform_01();
} while (x == 0 || x == 1);
return x;
#else
/* Rounding down. */
uint32_t z, x, hi, lo;
double s;
/*
* Draw an exponent, geometrically distributed, but give up if
* we get a run of more than 1088 zeros, which really means the
* system is broken.
*/
z = 0;
while ((x = crypto_rand_uint32()) == 0) {
if (z >= 1088)
/* Your bit sampler is broken. Go home. */
return __DBL_DENORM_MIN__;
z += 32;
}
z += clz32(x);
/*
* Pick the high and low chunks of a uniformly distributed
* significand in [2^52, 2^53).
*/
hi = (crypto_rand_uint32() & UINT32_C(0xfffff)) | 0x100000;
lo = crypto_rand_uint32();
/*
* Convert to a double in [2^52, 2^53). This is computed
* exactly, with no rounding, so 2^53 is impossible.
*/
s = hi*(double)4294967296 + lo;
/* Scale it into [1/2, 1) and then apply the exponent z. */
s = s * ldexp(1, -(53 + z));
/*
* Out of paranoia, clamp at the smallest subnormal. This
* should happen only with negligible probability, 2^-1075, so
* while it biases the distribution, you really have bigger
* problems than a negligibly biased distribution if this
* happens and I hope you bought a lot of lottery tickets.
* (Could do this earlier at the choice of z but I ran out of
* energy to figure out exactly the right fenceposts there.)
*/
return (s < __DBL_DENORM_MIN__ ? __DBL_DENORM_MIN__ : s);
#endif
}