Ticket #23061: uniform01.c

File uniform01.c, 10.2 KB (added by riastradh, 11 months ago)

uniform [0,1] sampler

Line 
1/*-
2 * Copyright (c) 2018 Taylor R. Campbell
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 * 1. Redistributions of source code must retain the above copyright
9 *    notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 *    notice, this list of conditions and the following disclaimer in the
12 *    documentation and/or other materials provided with the distribution.
13 *
14 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
15 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
17 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
18 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
20 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
21 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
22 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
23 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
24 * SUCH DAMAGE.
25 */
26
27/*
28 * If you take the uniform distribution on [0,1] or (0,1), i.e. the
29 * Lebesgue measure, and round it to the nearest floating-point number,
30 * the distribution on the result is the natural notion of the uniform
31 * distribution on the floating-point numbers in [0,1].
32 *
33 * Note that even if you round a sample from the open interval (0,1) of
34 * real numbers, the result might be rounded to 0 or 1, if it is in (0,
35 * 2^{emin - p - 1}] or [1 - eps/4, 1):
36 *
37 * - The interval near 0 is of negligible size in any format larger
38 *   than binary16 -- less than 2^-126 in binary32 even if you use
39 *   flush-to-zero, less than 2^-1022 in binary64 -- so you don't need
40 *   to worry about it and can safely compute, e.g., 1/x or log(x).
41 *
42 * - However, the probability of getting a real number in the interval
43 *   near 1 is 2^-54 in binary64, which is small but not negligible, so
44 *   I recommend not rejecting it unless you must avoid 1, for instance
45 *   because you need to compute log1p(-x).
46 *
47 * In infinite-precision floating-point, the exponent is geometrically
48 * distributed in {-1, -2, -3, ...}, and the fractional part of the
49 * significand is a uniformly distributed infinite bit string.  If we
50 * round this to finite precision, the exponent remains the same
51 * (unless rounded to subnormal or zero) and the significand is rounded
52 * to p bits.  The outcome of rounding is not actually affected by more
53 * than a finite number of significand bits, with one tricky exception,
54 * so you don't need to sample an infinite stream of bits and round it;
55 * to get this, you can:
56 *
57 * 1. flip a coin until you get one heads, to pick a geometrically
58 *    distributed exponent;
59 *
60 * 2. flip a coin >p+1 times to choose the uniformly distributed
61 *    fractional part of a significand, and round it to p bits with a
62 *    trick to avoid a bias toward `even'.
63 *
64 * If you drew p bits for the fractional part of the significand, there
65 * would be no rounding, and you would exclude from the possible
66 * significands the interval (2 - eps/2, 2] of real numbers that are
67 * rounded to 2 -- small but not nonnegligible.  If you draw >p bits
68 * and round it, then no real number in [1, 2] is excluded.
69 *
70 * However, there is a slight bias toward `even' outcomes -- that is,
71 * outcomes whose least significant bit is zero -- because there's a
72 * nonzero probability that the finite significand you choose is a tie,
73 * exactly between two floating-point numbers, which is broken by
74 * choosing the even option.  This event has probability zero in
75 * Lebesgue measure.  You can simply break the ties by drawing >p+1
76 * bits and forcing the least significant bit to be 1.
77 *
78 * The attached file uniform01.c implements this algorithm in the
79 * function sample_uniform_01, with parameters chosen for IEEE 754
80 * binary64 (`double-precision') floating-point.
81 *
82 * As it happens, the variant algorithm that chooses exactly p bits and
83 * _doesn't_ round them gives a uniform distribution on [0, 1 - eps/4),
84 * which incidentally effectively never returns zero, and so yields
85 * what you might think of `the uniform distribution on floating-point
86 * numbers in (0, 1)'.  I named it sample_uniform_01_open, and added
87 * some paranoia to clamp the result at the smallest subnormal in case
88 * your PRNG is broken.  But it also works to just do rejection
89 * sampling on sample_uniform_01 and reduce the amount of code to worry
90 * about.
91 *
92 * Note: This code does not run in constant time -- the time it takes
93 * is proportional to the magnitude of the exponent.  (In principle, it
94 * may also take variable time for subnormal arithmetic, but the
95 * probability of that is negligible for any format larger than
96 * binary16.)
97 *
98 * You could make it run in constant time by always flipping the
99 * maximum number of coins for an exponent, and counting the leading
100 * zero bits.  To improve performance, you can safely truncate the
101 * exponent at, say, 128 -- if you flip 128 coins and you get tails
102 * every time, either your coin flipper is broken or there's a glitch
103 * in the matrix.
104 */
105
106#include <float.h>
107#include <math.h>
108#include <stdint.h>
109
110/**
111 * Draw an integer in [0, 2^32) uniformly at random.
112 */
113#include <stdlib.h>             /* arc4random */
114uint32_t
115crypto_rand_uint32(void)
116{
117        return arc4random();    /* XXX */
118}
119
120/**
121 * Count number of one bits in 32-bit word
122 */
123unsigned
124bitcount32(uint32_t x)
125{
126
127        /* Count two-bit groups.  */
128        x -= (x >> 1) & UINT32_C(0x55555555);
129
130        /* Count four-bit groups.  */
131        x = ((x >> 2) & UINT32_C(0x33333333)) + (x & UINT32_C(0x33333333));
132
133        /* Count eight-bit groups.  */
134        x = (x + (x >> 4)) & UINT32_C(0x0f0f0f0f);
135
136        /* Sum all eight-bit groups, and extract the sum.  */
137        return (x * UINT32_C(0x01010101)) >> 24;
138}
139
140/**
141 * Count leading zeros in 32-bit word
142 */
143unsigned
144clz32(uint32_t x)
145{
146
147        /* Round up to a power of two.  */
148        x |= x >> 1;
149        x |= x >> 2;
150        x |= x >> 4;
151        x |= x >> 8;
152        x |= x >> 16;
153
154        /* Subtract count of one bits from 32.  */
155        return (32 - bitcount32(x));
156}
157
158/*
159 * The following sample_uniform_01 is tailored for IEEE 754 binary64
160 * floating-point or smaller.  It can be adapted to larger
161 * floating-point formats like i387 80-bit or IEEE 754 binary128, but
162 * it may require sampling more bits.
163 */
164CTASSERT(FLT_RADIX == 2);
165CTASSERT(-DBL_MIN_EXP <= 1021);
166CTASSERT(DBL_MANT_DIG <= 53);
167
168/**
169 * Sample a floating-point number in [0, 1] with uniform distribution.
170 *
171 * Note that the probability of returning 0 is less than 2^-1074, so
172 * callers need not check for it.  However, callers that cannot handle
173 * rounding to 1 must deal with that, because it occurs with
174 * probability 2^-54, which is small but nonnegligible.
175 */
176double
177sample_uniform_01(void)
178{
179        uint32_t z, x, hi, lo;
180        double s;
181
182        /*
183         * Draw an exponent, geometrically distributed, but give up if
184         * we get a run of more than 1088 zeros, which really means the
185         * system is broken.
186         */
187        z = 0;
188        while ((x = crypto_rand_uint32()) == 0) {
189                if (z >= 1088)
190                        /* Your bit sampler is broken.  Go home.  */
191                        return 0;
192                z += 32;
193        }
194        z += clz32(x);
195
196        /*
197         * Pick 32-bit halves of an odd normalized significand.
198         * Picking it odd breaks ties in the subsequent rounding, which
199         * occur only with measure zero in the uniform distribution on
200         * [0, 1].
201         */
202        hi = crypto_rand_uint32() | UINT32_C(0x80000000);
203        lo = crypto_rand_uint32() | UINT32_C(0x00000001);
204
205        /* Round to nearest scaled significand in [2^63, 2^64].  */
206        s = hi*(double)4294967296 + lo;
207
208        /* Rescale into [1/2, 1] and apply exponent in one swell foop.  */
209        return s * ldexp(1, -(64 + z));
210}
211
212/*
213 * The following sample_uniform_01_open is tailored for IEEE 754
214 * binary64 floating-point exactly.  It can be adapted to smaller or
215 * larger floating-point formats, but it must fit the bits exactly to
216 * avoid rounding.
217 */
218CTASSERT(FLT_RADIX == 2);
219CTASSERT(-DBL_MIN_EXP == 1021);
220CTASSERT(DBL_MANT_DIG == 53);
221
222/**
223 * Sample a floating-point number in (2^{emin-p-1}, 1 - eps/4) with
224 * uniform distribution.  The outcome lies in [2^{emin-p}, 1 - eps/2].
225 * This is sometimes described as `(0, 1)', but you are guaranteed not
226 * to get the floating-point rounding of any real number in the
227 * intervals [0, 2^{emin-p-1}) or (1 - eps/4, 1].
228 *
229 * Use this only if you are sure you can't handle the value 1, because
230 * the width of the interval that is rounded to 1, 2^-54 in binary64,
231 * is small but not negligible.  sample_uniform_01() is already
232 * effectively guaranteed not to return 0.
233 */
234double
235sample_uniform_01_open(void)
236{
237#if 0
238        /* Rejection sampling.  */
239        double x;
240
241        /* 0 is effectively impossible -- probability <2^-1074.  */
242        do {
243                x = uniform_01();
244        } while (x == 0 || x == 1);
245
246        return x;
247#else
248        /* Rounding down.  */
249        uint32_t z, x, hi, lo;
250        double s;
251
252        /*
253         * Draw an exponent, geometrically distributed, but give up if
254         * we get a run of more than 1088 zeros, which really means the
255         * system is broken.
256         */
257        z = 0;
258        while ((x = crypto_rand_uint32()) == 0) {
259                if (z >= 1088)
260                        /* Your bit sampler is broken.  Go home.  */
261                        return __DBL_DENORM_MIN__;
262                z += 32;
263        }
264        z += clz32(x);
265
266        /*
267         * Pick the high and low chunks of a uniformly distributed
268         * significand in [2^52, 2^53).
269         */
270        hi = (crypto_rand_uint32() & UINT32_C(0xfffff)) | 0x100000;
271        lo = crypto_rand_uint32();
272
273        /*
274         * Convert to a double in [2^52, 2^53).  This is computed
275         * exactly, with no rounding, so 2^53 is impossible.
276         */
277        s = hi*(double)4294967296 + lo;
278
279        /* Scale it into [1/2, 1) and then apply the exponent z.  */
280        s = s * ldexp(1, -(53 + z));
281
282        /*
283         * Out of paranoia, clamp at the smallest subnormal.  This
284         * should happen only with negligible probability, 2^-1075, so
285         * while it biases the distribution, you really have bigger
286         * problems than a negligibly biased distribution if this
287         * happens and I hope you bought a lot of lottery tickets.
288         * (Could do this earlier at the choice of z but I ran out of
289         * energy to figure out exactly the right fenceposts there.)
290         */
291        return (s < __DBL_DENORM_MIN__ ? __DBL_DENORM_MIN__ : s);
292#endif
293}