t/rng.c

   1 /*
   2  * Simple random number generator for testing.
   3  * Copyright © 2022-2024 Nick Bowler
   4  *
   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.
   9  *
  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.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program.  If not, see <https://www.gnu.org/licenses/>.
  17  *
  18  * The RNG implementation is adapted from xoshiro256+ by David Blackman
  19  * and Sebastiano Vigna, originally distributed under the Creative Commons
  20  * Zero public domain dedication.
  21  */
  22
  23 #include <config.h>
  24 #include <stdio.h>
  25 #include <stdlib.h>
  26 #include <string.h>
  27 #include <errno.h>
  28 #include <limits.h>
  29 #include <inttypes.h>
  30
  31 #if !TEST_RNG_NO_EXTERNAL_API
  32 #  include "test.h"
  33 #endif
  34
  35 #define B64(x) ((x) & 0xffffffffffffffffull)
  36 #define B32(x) ((x) & 0xffffffffu)
  37
  38 struct test_rng {
  39         unsigned long long state[4];
  40 };
  41
  42 /* Rotate val left by n bits.  The behaviour is undefined if n is zero. */
  43 static unsigned long long rot_left64(unsigned long long val, int n)
  44 {
  45         return B64( (val << n) | (val >> (64 - n)) );
  46 }
  47
  48 static unsigned long long xoshiro256p(unsigned long long *s)
  49 {
  50         unsigned long long tmp, ret;
  51
  52         ret = B64(s[0] + s[3]);
  53         tmp = B64(s[1] << 17);
  54
  55         s[2] ^= s[0];
  56         s[3] ^= s[1];
  57         s[1] ^= s[2];
  58         s[0] ^= s[3];
  59         s[2] ^= tmp;
  60         s[3] = rot_left64(s[3], 45);
  61
  62         return ret;
  63 }
  64
  65 /*
  66  * 32-bit random number generator for seed expansion.
  67  *
  68  * The output filter is from the public domain MurmurHash3 by Austin Appleby.
  69  *
  70  * The state is incremented by the fractional part of the golden ratio,
  71  * i.e., floor( (sqrt(5)-1)/2 * 2**32 )
  72  */
  73 static uint_least32_t splitmix32(uint_least32_t *state)
  74 {
  75         uint_least32_t z;
  76
  77         z = B32(*state += 0x9e3779b9);
  78         z = B32((z ^ (z >> 16)) * 0x85ebca6b);
  79         z = B32((z ^ (z >> 13)) * 0xc2b2ae35);
  80
  81         return z ^ (z >> 16);
  82 }
  83
  84 static uint_least64_t seed64(uint_least32_t *state)
  85 {
  86         uint_least64_t x;
  87
  88         x = splitmix32(state);
  89         return (x << 32) | splitmix32(state);
  90 }
  91
  92 #if !TEST_RNG_NO_EXTERNAL_API
  93 struct test_rng *test_rng_alloc(const char *seed_str)
  94 {
  95         test_uintmax seed_val;
  96         struct test_rng *rng;
  97         uint_least32_t seed;
  98
  99         if (!test_strtoumax(&seed_val, seed_str, 0xffffffff)) {
 100                 print_error("%s: invalid seed", seed_str);
 101                 return NULL;
 102         }
 103         seed = seed_val;
 104
 105         rng = malloc_nofail(sizeof *rng);
 106         rng->state[0] = seed64(&seed);
 107         rng->state[1] = seed64(&seed);
 108         rng->state[2] = seed64(&seed);
 109         rng->state[3] = seed64(&seed);
 110
 111         return rng;
 112 }
 113
 114 void test_rng_free(struct test_rng *rng)
 115 {
 116         free(rng);
 117 }
 118
 119 /* Calculate the least power of two greater than val, minus 1. */
 120 static unsigned rng_mask(unsigned val)
 121 {
 122         val |= val >> 1;
 123         val |= val >> 2;
 124         val |= val >> 4;
 125         val |= val >> 8;
 126
 127         if (UINT_MAX >= 65536)
 128                 val |= val >> 16;
 129
 130         return val;
 131 }
 132
 133 unsigned test_rng_uniform_int(struct test_rng *rng, unsigned max)
 134 {
 135         unsigned mask = rng_mask(max-1);
 136         unsigned long long val;
 137
 138         do {
 139                 val = ( xoshiro256p(rng->state) >> 32 ) & mask;
 140         } while (val >= max);
 141
 142         return val;
 143 }
 144 #endif