t/rng.c

   1 /*
   2  * Simple random number generator for testing.
   3  * Copyright © 2022-2023 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 #include "test.h"
  32
  33 #define B64(x) ((x) & 0xffffffffffffffff)
  34
  35 struct test_rng {
  36         unsigned long long state[4];
  37 };
  38
  39 /* Rotate val left by n bits.  The behaviour is undefined if n is zero. */
  40 static unsigned long long rot_left64(unsigned long long val, int n)
  41 {
  42         return B64( (val << n) | (val >> (64 - n)) );
  43 }
  44
  45 static unsigned long long xoshiro256p(unsigned long long *s)
  46 {
  47         unsigned long long tmp, ret;
  48
  49         ret = B64(s[0] + s[3]);
  50         tmp = B64(s[1] << 17);
  51
  52         s[2] ^= s[0];
  53         s[3] ^= s[1];
  54         s[1] ^= s[2];
  55         s[0] ^= s[3];
  56         s[2] ^= tmp;
  57         s[3] = rot_left64(s[3], 45);
  58
  59         return ret;
  60 }
  61
  62 static unsigned long long splitmix64(unsigned long long *state)
  63 {
  64         unsigned long long z;
  65
  66         z = B64(*state += 0x9e3779b97f4a7c15);
  67         z = B64((z ^ (z >> 30)) * 0xbf58476d1ce4e5b9);
  68         z = B64((z ^ (z >> 27)) * 0x94d049bb133111eb);
  69
  70         return z ^ (z >> 31);
  71 }
  72
  73 struct test_rng *test_rng_alloc(const char *seed_str)
  74 {
  75         unsigned long long seed;
  76         uintmax_t limit, seed_val;
  77         struct test_rng *rng;
  78
  79         limit  = (uintmax_t)0xffffffff;
  80         limit |= (limit << 16 << 16);
  81
  82         if (!test_strtoumax(&seed_val, seed_str, limit)) {
  83                 fprintf(stderr, "%s: invalid seed\n", seed_str);
  84                 return NULL;
  85         }
  86
  87         rng = malloc_nofail(sizeof *rng);
  88         rng->state[0] = splitmix64(&seed);
  89         rng->state[1] = splitmix64(&seed);
  90         rng->state[2] = splitmix64(&seed);
  91         rng->state[3] = splitmix64(&seed);
  92
  93         return rng;
  94 }
  95
  96 void test_rng_free(struct test_rng *rng)
  97 {
  98         free(rng);
  99 }
 100
 101 /* Calculate the least power of two greater than val, minus 1. */
 102 static unsigned rng_mask(unsigned val)
 103 {
 104         val |= val >> 1;
 105         val |= val >> 2;
 106         val |= val >> 4;
 107         val |= val >> 8;
 108
 109         if (UINT_MAX >= 65536)
 110                 val |= val >> 16;
 111
 112         return val;
 113 }
 114
 115 unsigned test_rng_uniform_int(struct test_rng *rng, unsigned max)
 116 {
 117         unsigned mask = rng_mask(max-1);
 118         unsigned long long val;
 119
 120         do {
 121                 val = ( xoshiro256p(rng->state) >> 32 ) & mask;
 122         } while (val >= max);
 123
 124         return val;
 125 }