1 /* $NetBSD: random.c,v 1.8 2025/01/26 16:25:38 christos Exp $ */ 2 3 /* 4 * Copyright (C) Internet Systems Consortium, Inc. ("ISC") 5 * 6 * SPDX-License-Identifier: MPL-2.0 7 * 8 * This Source Code Form is subject to the terms of the Mozilla Public 9 * License, v. 2.0. If a copy of the MPL was not distributed with this 10 * file, you can obtain one at https://mozilla.org/MPL/2.0/. 11 * 12 * See the COPYRIGHT file distributed with this work for additional 13 * information regarding copyright ownership. 14 */ 15 16 /* 17 * Portions of isc_random_uniform(): 18 * 19 * Copyright (c) 1996, David Mazieres <dm@uun.org> 20 * Copyright (c) 2008, Damien Miller <djm@openbsd.org> 21 * 22 * Permission to use, copy, modify, and distribute this software for any 23 * purpose with or without fee is hereby granted, provided that the above 24 * copyright notice and this permission notice appear in all copies. 25 * 26 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 27 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 28 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 29 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 30 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 31 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 32 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 33 */ 34 35 #include <inttypes.h> 36 #include <stdlib.h> 37 #include <string.h> 38 #include <unistd.h> 39 40 #include <isc/entropy.h> 41 #include <isc/random.h> 42 #include <isc/result.h> 43 #include <isc/thread.h> 44 #include <isc/types.h> 45 #include <isc/util.h> 46 47 /* 48 * Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org) 49 * 50 * To the extent possible under law, the author has dedicated all 51 * copyright and related and neighboring rights to this software to the 52 * public domain worldwide. This software is distributed without any 53 * warranty. 54 * 55 * See <http://creativecommons.org/publicdomain/zero/1.0/>. 56 */ 57 58 /* 59 * This is xoshiro128** 1.0, our 32-bit all-purpose, rock-solid generator. 60 * It has excellent (sub-ns) speed, a state size (128 bits) that is large 61 * enough for mild parallelism, and it passes all tests we are aware of. 62 * 63 * The state must be seeded so that it is not everywhere zero. 64 */ 65 66 static thread_local bool initialized = false; 67 static thread_local uint32_t seed[4] = { 0 }; 68 69 static uint32_t 70 rotl(const uint32_t x, int k) { 71 return (x << k) | (x >> (32 - k)); 72 } 73 74 static uint32_t 75 next(void) { 76 uint32_t result_starstar, t; 77 78 result_starstar = rotl(seed[0] * 5, 7) * 9; 79 t = seed[1] << 9; 80 81 seed[2] ^= seed[0]; 82 seed[3] ^= seed[1]; 83 seed[1] ^= seed[2]; 84 seed[0] ^= seed[3]; 85 86 seed[2] ^= t; 87 88 seed[3] = rotl(seed[3], 11); 89 90 return result_starstar; 91 } 92 93 static void 94 isc__random_initialize(void) { 95 if (initialized) { 96 return; 97 } 98 99 #if FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION 100 /* 101 * A fixed seed helps with problem reproduction when fuzzing. It must be 102 * non-zero else xoshiro128starstar will generate only zeroes, and the 103 * first result needs to be non-zero as expected by random_test.c 104 */ 105 seed[0] = 1; 106 #endif /* if FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION */ 107 108 while (seed[0] == 0 && seed[1] == 0 && seed[2] == 0 && seed[3] == 0) { 109 isc_entropy_get(seed, sizeof(seed)); 110 } 111 initialized = true; 112 } 113 114 uint8_t 115 isc_random8(void) { 116 isc__random_initialize(); 117 return (uint8_t)next(); 118 } 119 120 uint16_t 121 isc_random16(void) { 122 isc__random_initialize(); 123 return (uint16_t)next(); 124 } 125 126 uint32_t 127 isc_random32(void) { 128 isc__random_initialize(); 129 return next(); 130 } 131 132 void 133 isc_random_buf(void *buf, size_t buflen) { 134 REQUIRE(buf != NULL); 135 REQUIRE(buflen > 0); 136 137 int i; 138 uint32_t r; 139 140 isc__random_initialize(); 141 142 for (i = 0; i + sizeof(r) <= buflen; i += sizeof(r)) { 143 r = next(); 144 memmove((uint8_t *)buf + i, &r, sizeof(r)); 145 } 146 r = next(); 147 memmove((uint8_t *)buf + i, &r, buflen % sizeof(r)); 148 return; 149 } 150 151 uint32_t 152 isc_random_uniform(uint32_t limit) { 153 isc__random_initialize(); 154 155 /* 156 * Daniel Lemire's nearly-divisionless unbiased bounded random numbers. 157 * 158 * https://lemire.me/blog/?p=17551 159 * 160 * The raw random number generator `next()` returns a 32-bit value. 161 * We do a 64-bit multiply `next() * limit` and treat the product as a 162 * 32.32 fixed-point value less than the limit. Our result will be the 163 * integer part (upper 32 bits), and we will use the fraction part 164 * (lower 32 bits) to determine whether or not we need to resample. 165 */ 166 uint64_t num = (uint64_t)next() * (uint64_t)limit; 167 /* 168 * In the fast path, we avoid doing a division in most cases by 169 * comparing the fraction part of `num` with the limit, which is 170 * a slight over-estimate for the exact resample threshold. 171 */ 172 if ((uint32_t)(num) < limit) { 173 /* 174 * We are in the slow path where we re-do the approximate test 175 * more accurately. The exact threshold for the resample loop 176 * is the remainder after dividing the raw RNG limit `1 << 32` 177 * by the caller's limit. We use a trick to calculate it 178 * within 32 bits: 179 * 180 * (1 << 32) % limit 181 * == ((1 << 32) - limit) % limit 182 * == (uint32_t)(-limit) % limit 183 * 184 * This division is safe: we know that `limit` is strictly 185 * greater than zero because of the slow-path test above. 186 */ 187 uint32_t residue = (uint32_t)(-limit) % limit; 188 /* 189 * Unless we get one of `N = (1 << 32) - residue` valid 190 * values, we reject the sample. This `N` is a multiple of 191 * `limit`, so our results will be unbiased; and `N` is the 192 * largest multiple that fits in 32 bits, so rejections are as 193 * rare as possible. 194 * 195 * There are `limit` possible values for the integer part of 196 * our fixed-point number. Each one corresponds to `N/limit` 197 * or `N/limit + 1` possible fraction parts. For our result to 198 * be unbiased, every possible integer part must have the same 199 * number of possible valid fraction parts. So, when we get 200 * the superfluous value in the `N/limit + 1` cases, we need 201 * to reject and resample. 202 * 203 * Because of the multiplication, the possible values in the 204 * fraction part are equally spaced by `limit`, with varying 205 * gaps at each end of the fraction's 32-bit range. We will 206 * choose a range of size `N` (a multiple of `limit`) into 207 * which valid fraction values must fall, with the rest of the 208 * 32-bit range covered by the `residue`. Lemire's paper says 209 * that exactly `N/limit` possible values spaced apart by 210 * `limit` will fit into our size `N` valid range, regardless 211 * of the size of the end gaps, the phase alignment of the 212 * values, or the position of the range. 213 * 214 * So, when a fraction value falls in the `residue` outside 215 * our valid range, it is superfluous, and we resample. 216 */ 217 while ((uint32_t)(num) < residue) { 218 num = (uint64_t)next() * (uint64_t)limit; 219 } 220 } 221 /* 222 * Return the integer part (upper 32 bits). 223 */ 224 return (uint32_t)(num >> 32); 225 } 226