2 * Simple random number generator for testing.
3 * Copyright © 2022 Nick Bowler
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 * This implementation is adapted from xoshiro256+ and splitmix64
19 * by David Blackman and Sebastiano Vigna, originally distributed
20 * under the Creative Commons Zero public domain dedication.
35 #define B64(x) ((x) & 0xffffffffffffffff)
38 unsigned long long state[4];
41 /* Rotate val left by n bits. The behaviour is undefined if n is zero. */
42 static unsigned long long rot_left64(unsigned long long val, int n)
44 return B64( (val << n) | (val >> (64 - n)) );
47 static unsigned long long xoshiro256p(unsigned long long *s)
49 unsigned long long ret = s[0];
59 s[3] = rot_left64(s[3], 45);
64 static unsigned long long splitmix64(unsigned long long *state)
68 z = B64(*state += 0x9e3779b97f4a7c15);
69 z = B64((z ^ (z >> 30)) * 0xbf58476d1ce4e5b9);
70 z = B64((z ^ (z >> 27)) * 0x94d049bb133111eb);
76 # define STRTOULL strtoull
78 /* HP-UX 11 has __strtoull in <inttypes.h> */
79 # define STRTOULL __strtoull
82 * Just fall back to strtoul -- in the worst case we just lose the ability
83 * to set all 64 bits of the seed.
85 # define STRTOULL strtoul
88 struct test_rng *test_rng_alloc(const char *seed_str)
90 unsigned long long seed;
95 seed = STRTOULL(seed_str, &end, 0);
97 fprintf(stderr, "%s: invalid seed\n", seed_str);
99 } else if (errno != 0) {
100 fprintf(stderr, "%s: invalid seed: %s\n",
101 seed_str, strerror(errno));
105 rng = malloc_nofail(sizeof *rng);
106 rng->state[0] = splitmix64(&seed);
107 rng->state[1] = splitmix64(&seed);
108 rng->state[2] = splitmix64(&seed);
109 rng->state[3] = splitmix64(&seed);
114 void test_rng_free(struct test_rng *rng)
119 #define BITS_PER_FP_DIGIT ( FLT_RADIX < 4 ? 1 \
120 : FLT_RADIX < 8 ? 2 \
121 : FLT_RADIX < 16 ? 3 \
122 : FLT_RADIX < 32 ? 4 \
125 double test_rng_uniform(struct test_rng *rng)
127 unsigned long long val = xoshiro256p(rng->state);
128 int prec = MIN(64, DBL_MANT_DIG*BITS_PER_FP_DIGIT);
130 return ldexp(val >> (64-prec), -prec);
133 /* Calculate the least power of two greater than val, minus 1. */
134 static unsigned rng_mask(unsigned val)
141 if (UINT_MAX >= 65536)
147 unsigned test_rng_uniform_int(struct test_rng *rng, unsigned max)
149 unsigned mask = rng_mask(max-1);
150 unsigned long long val;
153 val = ( xoshiro256p(rng->state) >> 32 ) & mask;
154 } while (val >= max);