t/xos256ss.c

   1 /*  Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org)
   2
   3 To the extent possible under law, the author has dedicated all copyright
   4 and related and neighboring rights to this software to the public domain
   5 worldwide. This software is distributed without any warranty.
   6
   7 See <http://creativecommons.org/publicdomain/zero/1.0/>. */
   8
   9 #include <stdint.h>
  10
  11 /* This is xoshiro256** 1.0, one of our all-purpose, rock-solid
  12    generators. It has excellent (sub-ns) speed, a state (256 bits) that is
  13    large enough for any parallel application, and it passes all tests we
  14    are aware of.
  15
  16    For generating just floating-point numbers, xoshiro256+ is even faster.
  17
  18    The state must be seeded so that it is not everywhere zero. If you have
  19    a 64-bit seed, we suggest to seed a splitmix64 generator and use its
  20    output to fill s. */
  21
  22 static inline uint64_t rotl(const uint64_t x, int k) {
  23         return (x << k) | (x >> (64 - k));
  24 }
  25
  26
  27 static uint64_t s[4];
  28
  29 uint64_t next(void) {
  30         const uint64_t result = rotl(s[1] * 5, 7) * 9;
  31
  32         const uint64_t t = s[1] << 17;
  33
  34         s[2] ^= s[0];
  35         s[3] ^= s[1];
  36         s[1] ^= s[2];
  37         s[0] ^= s[3];
  38
  39         s[2] ^= t;
  40
  41         s[3] = rotl(s[3], 45);
  42
  43         return result;
  44 }
  45
  46
  47 /* This is the jump function for the generator. It is equivalent
  48    to 2^128 calls to next(); it can be used to generate 2^128
  49    non-overlapping subsequences for parallel computations. */
  50
  51 void jump(void) {
  52         static const uint64_t JUMP[] = { 0x180ec6d33cfd0aba, 0xd5a61266f0c9392c, 0xa9582618e03fc9aa, 0x39abdc4529b1661c };
  53
  54         uint64_t s0 = 0;
  55         uint64_t s1 = 0;
  56         uint64_t s2 = 0;
  57         uint64_t s3 = 0;
  58         for(int i = 0; i < sizeof JUMP / sizeof *JUMP; i++)
  59                 for(int b = 0; b < 64; b++) {
  60                         if (JUMP[i] & UINT64_C(1) << b) {
  61                                 s0 ^= s[0];
  62                                 s1 ^= s[1];
  63                                 s2 ^= s[2];
  64                                 s3 ^= s[3];
  65                         }
  66                         next();
  67                 }
  68
  69         s[0] = s0;
  70         s[1] = s1;
  71         s[2] = s2;
  72         s[3] = s3;
  73 }
  74
  75
  76
  77 /* This is the long-jump function for the generator. It is equivalent to
  78    2^192 calls to next(); it can be used to generate 2^64 starting points,
  79    from each of which jump() will generate 2^64 non-overlapping
  80    subsequences for parallel distributed computations. */
  81
  82 void long_jump(void) {
  83         static const uint64_t LONG_JUMP[] = { 0x76e15d3efefdcbbf, 0xc5004e441c522fb3, 0x77710069854ee241, 0x39109bb02acbe635 };
  84
  85         uint64_t s0 = 0;
  86         uint64_t s1 = 0;
  87         uint64_t s2 = 0;
  88         uint64_t s3 = 0;
  89         for(int i = 0; i < sizeof LONG_JUMP / sizeof *LONG_JUMP; i++)
  90                 for(int b = 0; b < 64; b++) {
  91                         if (LONG_JUMP[i] & UINT64_C(1) << b) {
  92                                 s0 ^= s[0];
  93                                 s1 ^= s[1];
  94                                 s2 ^= s[2];
  95                                 s3 ^= s[3];
  96                         }
  97                         next();
  98                 }
  99
 100         s[0] = s0;
 101         s[1] = s1;
 102         s[2] = s2;
 103         s[3] = s3;
 104 }