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+ and splitmix64
  19  * by David Blackman and Sebastiano Vigna, originally distributed under
  20  * the Creative Commons Zero public domain dedication.
  21  */
  22
  23 #include <config.h>
  24 #include <stdio.h>
  25 #include <stdlib.h>
  26 #include <string.h>
  27 #include <inttypes.h>
  28 #include <errno.h>
  29 #include <limits.h>
  30
  31 #if !TEST_RNG_NO_EXTERNAL_API
  32 #  include "test.h"
  33 #endif
  34
  35 #define B64(x) ((x) & 0xffffffffffffffff)
  36
  37 struct test_rng {
  38         unsigned long long state[4];
  39 };
  40
  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)
  43 {
  44         return B64( (val << n) | (val >> (64 - n)) );
  45 }
  46
  47 static unsigned long long xoshiro256p(unsigned long long *s)
  48 {
  49         unsigned long long tmp, ret;
  50
  51         ret = B64(s[0] + s[3]);
  52         tmp = B64(s[1] << 17);
  53
  54         s[2] ^= s[0];
  55         s[3] ^= s[1];
  56         s[1] ^= s[2];
  57         s[0] ^= s[3];
  58         s[2] ^= tmp;
  59         s[3] = rot_left64(s[3], 45);
  60
  61         return ret;
  62 }
  63
  64 static unsigned long long splitmix64(unsigned long long *state)
  65 {
  66         unsigned long long z;
  67
  68         z = B64(*state += 0x9e3779b97f4a7c15);
  69         z = B64((z ^ (z >> 30)) * 0xbf58476d1ce4e5b9);
  70         z = B64((z ^ (z >> 27)) * 0x94d049bb133111eb);
  71
  72         return z ^ (z >> 31);
  73 }
  74
  75 #if !TEST_RNG_NO_EXTERNAL_API
  76 struct test_rng *test_rng_alloc(const char *seed_str)
  77 {
  78         unsigned long long seed;
  79         uintmax_t limit, seed_val;
  80         struct test_rng *rng;
  81
  82         limit  = (uintmax_t)0xffffffff;
  83         limit |= (limit << 16 << 16);
  84
  85         if (!test_strtoumax(&seed_val, seed_str, limit)) {
  86                 print_error("%s: invalid seed", seed_str);
  87                 return NULL;
  88         }
  89         seed = seed_val;
  90
  91         rng = malloc_nofail(sizeof *rng);
  92         rng->state[0] = splitmix64(&seed);
  93         rng->state[1] = splitmix64(&seed);
  94         rng->state[2] = splitmix64(&seed);
  95         rng->state[3] = splitmix64(&seed);
  96
  97         return rng;
  98 }
  99
 100 void test_rng_free(struct test_rng *rng)
 101 {
 102         free(rng);
 103 }
 104
 105 /* Calculate the least power of two greater than val, minus 1. */
 106 static unsigned rng_mask(unsigned val)
 107 {
 108         val |= val >> 1;
 109         val |= val >> 2;
 110         val |= val >> 4;
 111         val |= val >> 8;
 112
 113         if (UINT_MAX >= 65536)
 114                 val |= val >> 16;
 115
 116         return val;
 117 }
 118
 119 unsigned test_rng_uniform_int(struct test_rng *rng, unsigned max)
 120 {
 121         unsigned mask = rng_mask(max-1);
 122         unsigned long long val;
 123
 124         do {
 125                 val = ( xoshiro256p(rng->state) >> 32 ) & mask;
 126         } while (val >= max);
 127
 128         return val;
 129 }
 130 #endif