NanoMagMC  v0.2
Monte Carlo Simulation Software for Atomistic Models of Magnetic Materials
xoroshiro128.hpp
Go to the documentation of this file.
1 #pragma once
2 
3 //
4 // The code in this file is free and unencumbered software released
5 // into the public domain.
6 //
7 // <http://creativecommons.org/publicdomain/zero/1.0/>.
8 //
9 
10 #ifndef RNG_H
11 #define RNG_H
12 
13 #include <cstdint>
14 #include <random>
15 #include <thread>
16 
17 #if defined(_MSC_VER)
18 /* Microsoft C/C++-compatible compiler
19  * Adaptive implementations
20  */
21 #include <intrin.h>
22 #define _rdtsc() __rdtsc()
23 
24 #elif defined(__GNUC__) && (defined(__x86_64__) || defined(__i386__))
25 /* GCC-compatible compiler, targeting x86/x86-64
26  * NO implementation needed
27  */
28 #include <x86intrin.h>
29 
30 #elif defined(__GNUC__) && defined(__ARM_NEON__)
31 /* GCC-compatible compiler, targeting ARM with NEON
32  * Adaptive implementations
33  */
34 #include <arm_neon.h>
35 
36 #elif defined(__GNUC__) && defined(__IWMMXT__)
37 /* GCC-compatible compiler, targeting ARM with WMMX
38  * Adaptive implementations
39  */
40 #include <mmintrin.h>
41 
42 #elif (defined(__GNUC__) || defined(__xlC__)) && (defined(__VEC__) || defined(__ALTIVEC__))
43 /* XLC or GCC-compatible compiler, targeting PowerPC with VMX/VSX
44  * Adaptive implementations
45  */
46 #include <altivec.h>
47 
48 #elif defined(__GNUC__) && defined(__SPE__)
49 /* GCC-compatible compiler, targeting PowerPC with SPE
50  * Adaptive implementations
51  */
52 #include <spe.h>
53 
54 #endif
55 
56 namespace rng {
57 
58 //
59 // Bitwise circular left shift.
60 //
61 static inline std::uint64_t
62 rotl(const std::uint64_t x, int k)
63 {
64  return (x << k) | (x >> (64 - k));
65 }
66 
67 //
68 // A simple random seed generator based on the entropy coming from the
69 // system thread/process scheduler. This is rather slow but seeds are
70 // normally generated very infrequently.
71 //
72 struct tsc_seed
73 {
74  using result_type = std::uint64_t;
75 
77  {
78  std::uint64_t base = _rdtsc();
79  std::uint64_t seed = base & 0xff;
80  for (int i = 1; i < 8; i++) {
81  std::this_thread::yield();
82  seed |= ((_rdtsc() - base) & 0xff) << (i << 3);
83  }
84  return seed;
85  }
86 };
87 
88 //
89 // A random seed generator based on std::random_device.
90 //
92 {
93  using result_type = std::uint64_t;
94 
96  {
97  std::random_device rd;
98  if (sizeof(result_type) > sizeof(std::random_device::result_type))
99  return rd() | (result_type{rd()} << 32);
100  else
101  return rd();
102  }
103 };
104 
105 //
106 // A random number generator with 64-bit internal state. It is based on
107 // the code from here: http://xoroshiro.di.unimi.it/splitmix64.c
108 //
109 struct rng64
110 {
111  using result_type = std::uint64_t;
112 
113  std::uint64_t state;
114 
115  rng64(std::uint64_t seed = 1) : state{seed} {}
116 
118  {
119  std::uint64_t z = (state += UINT64_C(0x9E3779B97F4A7C15));
120  z = (z ^ (z >> 30)) * UINT64_C(0xBF58476D1CE4E5B9);
121  z = (z ^ (z >> 27)) * UINT64_C(0x94D049BB133111EB);
122  return z ^ (z >> 31);
123  }
124 };
125 
126 //
127 // A random number generator with 128-bit internal state. It is based on
128 // the code from here: http://xoroshiro.di.unimi.it/xoroshiro128plus.c
129 //
130 struct rng128
131 {
132  using result_type = std::uint64_t;
133 
134  std::uint64_t state[2];
135 
136  rng128(std::uint64_t seed[2]) : state{seed[0], seed[1]} {}
137 
138  rng128(std::uint64_t s0, std::uint64_t s1) : state{s0, s1} {}
139 
140  rng128(std::uint64_t seed = 1)
141  {
142  rng64 seeder(seed);
143  state[0] = seeder();
144  state[1] = seeder();
145  }
146 
147  static constexpr result_type min() { return 0; }
148  static constexpr result_type max() { return ~ result_type(0); }
149 
151  {
152  const uint64_t s0 = state[0];
153  uint64_t s1 = state[1];
154  const uint64_t value = s0 + s1;
155 
156  s1 ^= s0;
157  state[0] = rotl(s0, 55) ^ s1 ^ (s1 << 14);
158  state[1] = rotl(s1, 36);
159 
160  return value;
161  }
162 
163  // This is the jump function for the generator. It is equivalent
164  // to 2 ^ 64 calls to next(); it can be used to generate 2 ^ 64
165  // non-overlapping subsequences for parallel computations.
166  void jump()
167  {
168  static const std::uint64_t j[] = {0xbeac0467eba5facb, 0xd86b048b86aa9922};
169 
170  std::uint64_t s0 = 0, s1 = 0;
171  for (std::size_t i = 0; i < sizeof j / sizeof j[0]; i++) {
172  for (int b = 0; b < 64; b++) {
173  if ((j[i] & UINT64_C(1) << b) != 0) {
174  s0 ^= state[0];
175  s1 ^= state[1];
176  }
177  operator()();
178  }
179  }
180 
181  state[0] = s0;
182  state[1] = s1;
183  }
184 };
185 
186 } // namespace rng
187 
188 #endif // RNG_H
result_type operator()()
Definition: xoroshiro128.hpp:150
result_type operator()()
Definition: xoroshiro128.hpp:95
rng128(std::uint64_t seed[2])
Definition: xoroshiro128.hpp:136
void jump()
Definition: xoroshiro128.hpp:166
Definition: xoroshiro128.hpp:109
Definition: xoroshiro128.hpp:130
result_type operator()()
Definition: xoroshiro128.hpp:117
result_type operator()()
Definition: xoroshiro128.hpp:76
static constexpr result_type max()
Definition: xoroshiro128.hpp:148
std::uint64_t result_type
Definition: xoroshiro128.hpp:93
Definition: xoroshiro128.hpp:91
std::uint64_t result_type
Definition: xoroshiro128.hpp:132
Definition: state.hpp:12
static constexpr result_type min()
Definition: xoroshiro128.hpp:147
std::uint64_t result_type
Definition: xoroshiro128.hpp:74
rng128(std::uint64_t seed=1)
Definition: xoroshiro128.hpp:140
std::uint64_t result_type
Definition: xoroshiro128.hpp:111
Definition: xoroshiro128.hpp:56
rng64(std::uint64_t seed=1)
Definition: xoroshiro128.hpp:115
rng128(std::uint64_t s0, std::uint64_t s1)
Definition: xoroshiro128.hpp:138
Definition: xoroshiro128.hpp:72
std::uint64_t state
Definition: xoroshiro128.hpp:113