test/rng.c

   1 /*
   2  *  Simple random number generator for testing.
   3  *  Copyright © 2022 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 <http://www.gnu.org/licenses/>.
  17  *
  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.
  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 <float.h>
  30 #include <math.h>
  31
  32 #include "test.h"
  33
  34 #define B64(x) ((x) & 0xffffffffffffffff)
  35
  36 struct test_rng {
  37         unsigned long long state[4];
  38 };
  39
  40 /* Rotate val left by n bits.  The behaviour is undefined if n is zero. */
  41 static unsigned long long rot_left64(unsigned long long val, int n)
  42 {
  43         return B64( (val << n) | (val >> (64 - n)) );
  44 }
  45
  46 static unsigned long long xoshiro256p(unsigned long long *s)
  47 {
  48         unsigned long long ret = s[0];
  49         unsigned long long t;
  50
  51         t = B64(s[1] << 17);
  52
  53         s[2] ^= s[0];
  54         s[3] ^= s[1];
  55         s[1] ^= s[2];
  56         s[0] ^= s[3];
  57         s[2] ^= t;
  58         s[3] = rot_left64(s[3], 45);
  59
  60         return ret;
  61 }
  62
  63 static unsigned long long splitmix64(unsigned long long *state)
  64 {
  65         unsigned long long z;
  66
  67         z = B64(*state += 0x9e3779b97f4a7c15);
  68         z = B64((z ^ (z >> 30)) * 0xbf58476d1ce4e5b9);
  69         z = B64((z ^ (z >> 27)) * 0x94d049bb133111eb);
  70
  71         return z ^ (z >> 31);
  72 }
  73
  74 struct test_rng *test_rng_alloc(const char *seed_str)
  75 {
  76         unsigned long long seed;
  77         struct test_rng *rng;
  78         char *end;
  79
  80         errno = 0;
  81         seed = strtoull(seed_str, &end, 0);
  82         if (*end != 0) {
  83                 fprintf(stderr, "%s: invalid seed\n", seed_str);
  84                 return NULL;
  85         } else if (errno != 0) {
  86                 fprintf(stderr, "%s: invalid seed: %s\n",
  87                                 seed_str, strerror(errno));
  88                 return NULL;
  89         }
  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 #define BITS_PER_FP_DIGIT ( FLT_RADIX <  4 ? 1 \
 106                           : FLT_RADIX <  8 ? 2 \
 107                           : FLT_RADIX < 16 ? 3 \
 108                           : FLT_RADIX < 32 ? 4 \
 109                           : 5 )
 110
 111 double test_rng_uniform(struct test_rng *rng)
 112 {
 113         unsigned long long val = xoshiro256p(rng->state);
 114         int prec = MIN(64, DBL_MANT_DIG*BITS_PER_FP_DIGIT);
 115
 116         return ldexp(val >> (64-prec), -prec);
 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 }