source/kit/mersenne_twister_64.c


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

#include "mersenne_twister_64.h"

#include "time.h"

void kit_mt64_init(kit_mt64_state_t *state, uint64_t seed) {
  state->mt[0] = seed;
  for (state->index = 1; state->index < KIT_MT64_N; state->index++)
    state->mt[state->index] = (6364136223846793005ull *
                                   (state->mt[state->index - 1] ^
                                    (state->mt[state->index - 1] >>
                                     62u)) +
                               state->index);
}

uint64_t kit_mt64_generate(kit_mt64_state_t *state) {
  static uint64_t const mag01[2] = { 0ull, 0xB5026F5AA96619E9ull };

  int      i;
  uint64_t x;

  if (state->index >= KIT_MT64_N) {
    for (i = 0; i < KIT_MT64_N - 156; i++) {
      x = (state->mt[i] & 0xFFFFFFFF80000000ull) |
          (state->mt[i + 1] & 0x7FFFFFFFull);
      state->mt[i] = state->mt[i + 156] ^ (x >> 1u) ^
                     mag01[(int) (x & 1ull)];
    }

    for (; i < KIT_MT64_N - 1; i++) {
      x = (state->mt[i] & 0xFFFFFFFF80000000ull) |
          (state->mt[i + 1] & 0x7FFFFFFFull);
      state->mt[i] = state->mt[i + (156 - KIT_MT64_N)] ^ (x >> 1u) ^
                     mag01[(int) (x & 1ull)];
    }

    x = (state->mt[KIT_MT64_N - 1] & 0xFFFFFFFF80000000ull) |
        (state->mt[0] & 0x7FFFFFFFull);
    state->mt[KIT_MT64_N - 1] = state->mt[156 - 1] ^ (x >> 1u) ^
                                mag01[(int) (x & 1ull)];

    state->index = 0;
  }

  x = state->mt[state->index++];

  x ^= (x >> 29u) & 0x5555555555555555ull;
  x ^= (x << 17u) & 0x71D67FFFEDA60000ull;
  x ^= (x << 37u) & 0xFFF7EEE000000000ull;
  x ^= (x >> 43u);

  return x;
}

uint64_t kit_mt64_seed() {
  struct timespec t;
  timespec_get(&t, TIME_UTC);

  kit_mt64_state_t s;
  kit_mt64_init(&s, (uint64_t) t.tv_nsec);

  return kit_mt64_generate(&s);
}