c_cpp/lib/klib/krng.h

   1 #ifndef KRNG_H
   2 #define KRNG_H
   3
   4 typedef struct {
   5         uint64_t s[2];
   6 } krng_t;
   7
   8 static inline uint64_t kr_splitmix64(uint64_t x)
   9 {
  10         uint64_t z = (x += 0x9E3779B97F4A7C15ULL);
  11         z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ULL;
  12         z = (z ^ (z >> 27)) * 0x94D049BB133111EBULL;
  13         return z ^ (z >> 31);
  14 }
  15
  16 static inline uint64_t kr_rand_r(krng_t *r)
  17 {
  18         const uint64_t s0 = r->s[0];
  19         uint64_t s1 = r->s[1];
  20         const uint64_t result = s0 + s1;
  21         s1 ^= s0;
  22         r->s[0] = (s0 << 55 | s0 >> 9) ^ s1 ^ (s1 << 14);
  23         r->s[1] = s0 << 36 | s0 >> 28;
  24         return result;
  25 }
  26
  27 static inline void kr_jump_r(krng_t *r)
  28 {
  29         static const uint64_t JUMP[] = { 0xbeac0467eba5facbULL, 0xd86b048b86aa9922ULL };
  30         uint64_t s0 = 0, s1 = 0;
  31         int i, b;
  32         for (i = 0; i < 2; ++i)
  33                 for (b = 0; b < 64; b++) {
  34                         if (JUMP[i] & 1ULL << b)
  35                                 s0 ^= r->s[0], s1 ^= r->s[1];
  36                         kr_rand_r(r);
  37                 }
  38         r->s[0] = s0, r->s[1] = s1;
  39 }
  40
  41 static inline void kr_srand_r(krng_t *r, uint64_t seed)
  42 {
  43         r->s[0] = kr_splitmix64(seed);
  44         r->s[1] = kr_splitmix64(r->s[0]);
  45 }
  46
  47 static inline double kr_drand_r(krng_t *r)
  48 {
  49         union { uint64_t i; double d; } u;
  50         u.i = 0x3FFULL << 52 | kr_rand_r(r) >> 12;
  51         return u.d - 1.0;
  52 }
  53
  54 #endif