source: libcfa/src/bits/random.hfa@ e62802f

ADT ast-experimental
Last change on this file since e62802f was c2dfa56a, checked in by Peter A. Buhr <pabuhr@…>, 3 years ago

remove SKULLDUGGERY dealing with CFA typedef numbering problem
  • Property mode set to 100644
File size: 9.8 KB
Line 
1//
2// Cforall Version 1.0.0 Copyright (C) 2022 University of Waterloo
3//
4// The contents of this file are covered under the licence agreement in the
5// file "LICENCE" distributed with Cforall.
6//
7// random.hfa --
8//
9// Author : Peter A. Buhr
10// Created On : Fri Jan 14 07:18:11 2022
11// Last Modified By : Peter A. Buhr
12// Last Modified On : Mon Mar 20 10:01:40 2023
13// Update Count : 180
14//
15
16#pragma once
17
18#include <stdint.h> // uintXX_t
19
20#define GLUE2( x, y ) x##y
21#define GLUE( x, y ) GLUE2( x, y )
22
23// Set default PRNG for architecture size.
24#ifdef __x86_64__ // 64-bit architecture
25 // 64-bit generators
26 //#define LEHMER64
27 //#define XORSHIFT_12_25_27
28 #define XOSHIRO256PP
29 //#define KISS_64
30
31 // 32-bit generators
32 //#define XORSHIFT_6_21_7
33 #define XOSHIRO128PP
34#else // 32-bit architecture
35 // 64-bit generators
36 //#define XORSHIFT_13_7_17
37 #define XOSHIRO256PP
38
39 // 32-bit generators
40 //#define XORSHIFT_6_21_7
41 #define XOSHIRO128PP
42#endif // __x86_64__
43
44// Define C/CFA PRNG name and random-state.
45
46#ifdef XOSHIRO256PP
47#define PRNG_NAME_64 xoshiro256pp
48#define PRNG_STATE_64_T GLUE(PRNG_NAME_64,_t)
49typedef struct { uint64_t s0, s1, s2, s3; } PRNG_STATE_64_T;
50#endif // XOSHIRO256PP
51
52#ifdef XOSHIRO128PP
53#define PRNG_NAME_32 xoshiro128pp
54#define PRNG_STATE_32_T GLUE(PRNG_NAME_32,_t)
55typedef struct { uint32_t s0, s1, s2, s3; } PRNG_STATE_32_T;
56#endif // XOSHIRO128PP
57
58#ifdef LEHMER64
59#define PRNG_NAME_64 lehmer64
60#define PRNG_STATE_64_T __uint128_t
61#endif // LEHMER64
62
63#ifdef WYHASH64
64#define PRNG_NAME_64 wyhash64
65#define PRNG_STATE_64_T uint64_t
66#endif // LEHMER64
67
68#ifdef XORSHIFT_13_7_17
69#define PRNG_NAME_64 xorshift_13_7_17
70#define PRNG_STATE_64_T uint64_t
71#endif // XORSHIFT_13_7_17
72
73#ifdef XORSHIFT_6_21_7
74#define PRNG_NAME_32 xorshift_6_21_7
75#define PRNG_STATE_32_T uint32_t
76#endif // XORSHIFT_6_21_7
77
78#ifdef XORSHIFT_12_25_27
79#define PRNG_NAME_64 xorshift_12_25_27
80#define PRNG_STATE_64_T uint64_t
81#endif // XORSHIFT_12_25_27
82
83#ifdef KISS_64
84#define PRNG_NAME_64 kiss_64
85#define PRNG_STATE_64_T GLUE(PRNG_NAME_64,_t)
86typedef struct { uint64_t z, w, jsr, jcong; } PRNG_STATE_64_T;
87#endif // KISS_^64
88
89#ifdef XORWOW
90#define PRNG_NAME_32 xorwow
91#define PRNG_STATE_32_T GLUE(PRNG_NAME_32,_t)
92typedef struct { uint32_t a, b, c, d, counter; } PRNG_STATE_32_T;
93#endif // XOSHIRO128PP
94
95#define PRNG_SET_SEED_64 GLUE(PRNG_NAME_64,_set_seed)
96#define PRNG_SET_SEED_32 GLUE(PRNG_NAME_32,_set_seed)
97
98
99// Default PRNG used by runtime.
100#ifdef __x86_64__ // 64-bit architecture
101#define PRNG_NAME PRNG_NAME_64
102#define PRNG_STATE_T PRNG_STATE_64_T
103#else // 32-bit architecture
104#define PRNG_NAME PRNG_NAME_32
105#define PRNG_STATE_T PRNG_STATE_32_T
106#endif // __x86_64__
107
108#define PRNG_SET_SEED GLUE(PRNG_NAME,_set_seed)
109
110
111// ALL PRNG ALGORITHMS ARE OPTIMIZED SO THAT THE PRNG LOGIC CAN HAPPEN IN PARALLEL WITH THE USE OF THE RESULT.
112// Specifically, the current random state is copied for returning, before computing the next value. As a consequence,
113// the set_seed routine primes the PRNG by calling it with the state so the seed is not return as the first random
114// value.
115
116
117#ifdef __cforall // don't include in C code (invoke.h)
118
119// https://prng.di.unimi.it/xoshiro256starstar.c
120//
121// This is xoshiro256++ 1.0, one of our all-purpose, rock-solid generators. It has excellent (sub-ns) speed, a state
122// (256 bits) that is large enough for any parallel application, and it passes all tests we are aware of.
123//
124// For generating just floating-point numbers, xoshiro256+ is even faster.
125//
126// The state must be seeded so that it is not everywhere zero. If you have a 64-bit seed, we suggest to seed a
127// splitmix64 generator and use its output to fill s.
128
129#ifndef XOSHIRO256PP
130typedef struct { uint64_t s0, s1, s2, s3; } xoshiro256pp_t;
131#endif // ! XOSHIRO256PP
132
133static inline uint64_t xoshiro256pp( xoshiro256pp_t & rs ) with(rs) {
134 inline uint64_t rotl( const uint64_t x, int k ) {
135 return (x << k) | (x >> (64 - k));
136 } // rotl
137
138 const uint64_t result = rotl( s0 + s3, 23 ) + s0;
139 const uint64_t t = s1 << 17;
140
141 s2 ^= s0;
142 s3 ^= s1;
143 s1 ^= s2;
144 s0 ^= s3;
145 s2 ^= t;
146 s3 = rotl( s3, 45 );
147 return result;
148} // xoshiro256pp
149
150static inline void xoshiro256pp_set_seed( xoshiro256pp_t & state, uint64_t seed ) {
151 state = (xoshiro256pp_t){ seed, seed, seed, seed };
152 xoshiro256pp( state ); // prime
153} // xoshiro256pp_set_seed
154
155// https://prng.di.unimi.it/xoshiro128plusplus.c
156//
157// This is xoshiro128++ 1.0, one of our 32-bit all-purpose, rock-solid generators. It has excellent speed, a state size
158// (128 bits) that is large enough for mild parallelism, and it passes all tests we are aware of.
159//
160// For generating just single-precision (i.e., 32-bit) floating-point numbers, xoshiro128+ is even faster.
161//
162// The state must be seeded so that it is not everywhere zero.
163
164#ifndef XOSHIRO128PP
165typedef struct { uint32_t s0, s1, s2, s3; } xoshiro128pp_t;
166#endif // ! XOSHIRO128PP
167
168static inline uint32_t xoshiro128pp( xoshiro128pp_t & rs ) with(rs) {
169 inline uint32_t rotl( const uint32_t x, int k ) {
170 return (x << k) | (x >> (32 - k));
171 } // rotl
172
173 const uint32_t result = rotl( s0 + s3, 7 ) + s0;
174 const uint32_t t = s1 << 9;
175
176 s2 ^= s0;
177 s3 ^= s1;
178 s1 ^= s2;
179 s0 ^= s3;
180 s2 ^= t;
181 s3 = rotl( s3, 11 );
182 return result;
183} // xoshiro128pp
184
185static inline void xoshiro128pp_set_seed( xoshiro128pp_t & state, uint32_t seed ) {
186 state = (xoshiro128pp_t){ seed, seed, seed, seed };
187 xoshiro128pp( state ); // prime
188} // xoshiro128pp_set_seed
189
190#ifdef __SIZEOF_INT128__
191 //--------------------------------------------------
192 static inline uint64_t lehmer64( __uint128_t & state ) {
193 __uint128_t ret = state;
194 state *= 0x_da94_2042_e4dd_58b5;
195 return ret >> 64;
196 } // lehmer64
197
198 static inline void lehmer64_set_seed( __uint128_t & state, uint64_t seed ) {
199 // The seed needs to be coprime with the 2^64 modulus to get the largest period, so no factors of 2 in the seed.
200 state = seed;
201 lehmer64( state ); // prime
202 } // lehmer64_set_seed
203
204 //--------------------------------------------------
205 static inline uint64_t wyhash64( uint64_t & state ) {
206 uint64_t ret = state;
207 state += 0x_60be_e2be_e120_fc15;
208 __uint128_t tmp;
209 tmp = (__uint128_t) ret * 0x_a3b1_9535_4a39_b70d;
210 uint64_t m1 = (tmp >> 64) ^ tmp;
211 tmp = (__uint128_t)m1 * 0x_1b03_7387_12fa_d5c9;
212 uint64_t m2 = (tmp >> 64) ^ tmp;
213 return m2;
214 } // wyhash64
215
216 static inline void wyhash64_set_seed( uint64_t & state, uint64_t seed ) {
217 state = seed;
218 wyhash64( state ); // prime
219 } // wyhash64_set_seed
220#endif // __SIZEOF_INT128__
221
222//--------------------------------------------------
223static inline uint64_t xorshift_13_7_17( uint64_t & state ) {
224 uint64_t ret = state;
225 state ^= state << 13;
226 state ^= state >> 7;
227 state ^= state << 17;
228 return ret;
229} // xorshift_13_7_17
230
231static inline void xorshift_13_7_17_set_seed( uint64_t & state, uint64_t seed ) {
232 state = seed;
233 xorshift_13_7_17( state ); // prime
234} // xorshift_13_7_17_set_seed
235
236//--------------------------------------------------
237// Marsaglia shift-XOR PRNG with thread-local state
238// Period is 4G-1
239// 0 is absorbing and must be avoided
240// Low-order bits are not particularly random
241static inline uint32_t xorshift_6_21_7( uint32_t & state ) {
242 uint32_t ret = state;
243 state ^= state << 6;
244 state ^= state >> 21;
245 state ^= state << 7;
246 return ret;
247} // xorshift_6_21_7
248
249static inline void xorshift_6_21_7_set_seed( uint32_t & state, uint32_t seed ) {
250 state = seed;
251 xorshift_6_21_7( state ); // prime
252} // xorshift_6_21_7_set_seed
253
254//--------------------------------------------------
255// The state must be seeded with a nonzero value.
256static inline uint64_t xorshift_12_25_27( uint64_t & state ) {
257 uint64_t ret = state;
258 state ^= state >> 12;
259 state ^= state << 25;
260 state ^= state >> 27;
261 return ret * 0x_2545_F491_4F6C_DD1D;
262} // xorshift_12_25_27
263
264static inline void xorshift_12_25_27_set_seed( uint64_t & state, uint64_t seed ) {
265 state = seed;
266 xorshift_12_25_27( state ); // prime
267} // xorshift_12_25_27_set_seed
268
269//--------------------------------------------------
270// The state must be seeded with a nonzero value.
271#ifndef KISS_64
272typedef struct { uint64_t z, w, jsr, jcong; } kiss_64_t;
273#endif // ! KISS_64
274
275static inline uint64_t kiss_64( kiss_64_t & rs ) with(rs) {
276 kiss_64_t ret = rs;
277 z = 36969 * (z & 65535) + (z >> 16);
278 w = 18000 * (w & 65535) + (w >> 16);
279 jsr ^= (jsr << 13);
280 jsr ^= (jsr >> 17);
281 jsr ^= (jsr << 5);
282 jcong = 69069 * jcong + 1234567;
283 return (((ret.z << 16) + ret.w) ^ ret.jcong) + ret.jsr;
284} // kiss_64
285
286static inline void kiss_64_set_seed( kiss_64_t & rs, uint64_t seed ) with(rs) {
287 z = 1; w = 1; jsr = 4; jcong = seed;
288 kiss_64( rs ); // prime
289} // kiss_64_set_seed
290
291//--------------------------------------------------
292// The state array must be initialized to non-zero in the first four words.
293#ifndef XORWOW
294typedef struct { uint32_t a, b, c, d, counter; } xorwow_t;
295#endif // ! XORWOW
296
297static inline uint32_t xorwow( xorwow_t & rs ) with(rs) {
298 // Algorithm "xorwow" from p. 5 of Marsaglia, "Xorshift RNGs".
299 uint32_t ret = a + counter;
300 uint32_t t = d;
301
302 uint32_t const s = a;
303 d = c;
304 c = b;
305 b = s;
306
307 t ^= t >> 2;
308 t ^= t << 1;
309 t ^= s ^ (s << 4);
310 a = t;
311 counter += 362437;
312 return ret;
313} // xorwow
314
315static inline void xorwow_set_seed( xorwow_t & rs, uint32_t seed ) {
316 rs = (xorwow_t){ seed, seed, seed, seed, 0 };
317 xorwow( rs ); // prime
318} // xorwow_set_seed
319
320//--------------------------------------------------
321// Used in __tls_rand_fwd
322#define M (1_l64u << 48_l64u)
323#define A (25214903917_l64u)
324#define AI (18446708753438544741_l64u)
325#define C (11_l64u)
326#define D (16_l64u)
327
328// Bi-directional LCG random-number generator
329static inline uint32_t LCGBI_fwd( uint64_t & rs ) {
330 rs = (A * rs + C) & (M - 1);
331 return rs >> D;
332} // LCGBI_fwd
333
334static inline uint32_t LCGBI_bck( uint64_t & rs ) {
335 unsigned int r = rs >> D;
336 rs = AI * (rs - C) & (M - 1);
337 return r;
338} // LCGBI_bck
339
340#undef M
341#undef A
342#undef AI
343#undef C
344#undef D
345
346#endif // __cforall
Note: See TracBrowser for help on using the repository browser.