1 /* $OpenBSD: random.c,v 1.22 2023/02/18 08:52:39 miod Exp $ */
2 /* $NetBSD: random.c,v 1.3 1995/04/22 07:44:05 cgd Exp $ */
3
4 /*
5 * Copyright (c) 1994
6 * The Regents of the University of California. All rights reserved.
7 *
8 * This code is derived from software contributed to Berkeley by
9 * Guy Harris at Network Appliance Corp.
10 *
11 * Redistribution and use in source and binary forms, with or without
12 * modification, are permitted provided that the following conditions
13 * are met:
14 * 1. Redistributions of source code must retain the above copyright
15 * notice, this list of conditions and the following disclaimer.
16 * 2. Redistributions in binary form must reproduce the above copyright
17 * notice, this list of conditions and the following disclaimer in the
18 * documentation and/or other materials provided with the distribution.
19 * 3. Neither the name of the University nor the names of its contributors
20 * may be used to endorse or promote products derived from this software
21 * without specific prior written permission.
22 *
23 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
24 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
27 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
28 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
29 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
32 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
33 * SUCH DAMAGE.
34 */
35
36 /*-
37 * Copyright (c) 2014 Taylor R. Campbell
38 * All rights reserved.
39 *
40 * Redistribution and use in source and binary forms, with or without
41 * modification, are permitted provided that the following conditions
42 * are met:
43 * 1. Redistributions of source code must retain the above copyright
44 * notice, this list of conditions and the following disclaimer.
45 * 2. Redistributions in binary form must reproduce the above copyright
46 * notice, this list of conditions and the following disclaimer in the
47 * documentation and/or other materials provided with the distribution.
48 *
49 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
50 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
51 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
52 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
53 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
54 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
55 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
56 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
57 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
58 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
59 * SUCH DAMAGE.
60 */
61
62 #include <err.h>
63 #include <errno.h>
64 #include <math.h>
65 #include <stdio.h>
66 #include <stdlib.h>
67 #include <unistd.h>
68
69 __dead void usage(void);
70
71 int
clz64(uint64_t x)72 clz64(uint64_t x)
73 {
74 static const uint64_t mask[] = {
75 0xffffffff00000000ULL, 0xffff0000, 0xff00, 0xf0, 0xc, 0x2,
76 };
77 static const int zeroes[] = {
78 32, 16, 8, 4, 2, 1,
79 };
80 int clz = 0;
81 int i;
82
83 if (x == 0)
84 return 64;
85
86 for (i = 0; i < 6; i++) {
87 if ((x & mask[i]) == 0)
88 clz += zeroes[i];
89 else
90 x >>= zeroes[i];
91 }
92
93 return clz;
94 }
95
96 uint64_t
random64(void)97 random64(void)
98 {
99 uint64_t r;
100
101 arc4random_buf(&r, sizeof(uint64_t));
102
103 return r;
104 }
105
106 /*
107 * random_real: Generate a stream of bits uniformly at random and
108 * interpret it as the fractional part of the binary expansion of a
109 * number in [0, 1], 0.00001010011111010100...; then round it.
110 */
111 double
random_real(void)112 random_real(void)
113 {
114 int exponent = -64;
115 uint64_t significand;
116 int shift;
117
118 /*
119 * Read zeros into the exponent until we hit a one; the rest
120 * will go into the significand.
121 */
122 while (__predict_false((significand = random64()) == 0)) {
123 exponent -= 64;
124 /*
125 * If the exponent falls below -1074 = emin + 1 - p,
126 * the exponent of the smallest subnormal, we are
127 * guaranteed the result will be rounded to zero. This
128 * case is so unlikely it will happen in realistic
129 * terms only if random64 is broken.
130 */
131 if (__predict_false(exponent < -1074))
132 return 0;
133 }
134
135 /*
136 * There is a 1 somewhere in significand, not necessarily in
137 * the most significant position. If there are leading zeros,
138 * shift them into the exponent and refill the less-significant
139 * bits of the significand. Can't predict one way or another
140 * whether there are leading zeros: there's a fifty-fifty
141 * chance, if random64 is uniformly distributed.
142 */
143 shift = clz64(significand);
144 if (shift != 0) {
145 exponent -= shift;
146 significand <<= shift;
147 significand |= (random64() >> (64 - shift));
148 }
149
150 /*
151 * Set the sticky bit, since there is almost surely another 1
152 * in the bit stream. Otherwise, we might round what looks
153 * like a tie to even when, almost surely, were we to look
154 * further in the bit stream, there would be a 1 breaking the
155 * tie.
156 */
157 significand |= 1;
158
159 /*
160 * Finally, convert to double (rounding) and scale by
161 * 2^exponent.
162 */
163 return ldexp((double)significand, exponent);
164 }
165
166 int
main(int argc,char * argv[])167 main(int argc, char *argv[])
168 {
169 double denom;
170 int ch, random_exit, selected, unbuffer_output;
171 char *ep;
172
173 if (pledge("stdio", NULL) == -1)
174 err(1, "pledge");
175
176 random_exit = unbuffer_output = 0;
177 while ((ch = getopt(argc, argv, "erh")) != -1)
178 switch (ch) {
179 case 'e':
180 random_exit = 1;
181 break;
182 case 'r':
183 unbuffer_output = 1;
184 break;
185 default:
186 case 'h':
187 usage();
188 }
189
190 argc -= optind;
191 argv += optind;
192
193 switch (argc) {
194 case 0:
195 denom = 2;
196 break;
197 case 1:
198 errno = 0;
199 denom = strtod(*argv, &ep);
200 if (errno == ERANGE)
201 err(1, "%s", *argv);
202 if (denom < 1 || *ep != '\0')
203 errx(1, "denominator is not valid.");
204 break;
205 default:
206 usage();
207 }
208
209 /* Return a random exit status between 0 and min(denom - 1, 255). */
210 if (random_exit) {
211 if (denom > 256)
212 denom = 256;
213
214 return arc4random_uniform(denom);
215 }
216
217 /*
218 * Act as a filter, randomly choosing lines of the standard input
219 * to write to the standard output.
220 */
221 if (unbuffer_output)
222 setvbuf(stdout, NULL, _IONBF, 0);
223
224 /*
225 * Select whether to print the first line. (Prime the pump.)
226 * We find a random number between 0 and 1 and, if it's < 1 / denom,
227 * we select the line.
228 */
229 selected = random_real() < 1 / denom;
230 while ((ch = getchar()) != EOF) {
231 int retch;
232
233 if (selected) {
234 errno = 0;
235 retch = putchar(ch);
236 if (retch == EOF && errno)
237 err(2, "putchar");
238 }
239 if (ch == '\n')
240 /* Now see if the next line is to be printed. */
241 selected = random_real() < 1 / denom;
242 }
243 if (ferror(stdin))
244 err(2, "stdin");
245 return 0;
246 }
247
248 void
usage(void)249 usage(void)
250 {
251
252 (void)fprintf(stderr, "usage: %s [-er] [denominator]\n", getprogname());
253 exit(1);
254 }
255