source: libcfa/src/bits/random.hfa@ 814a4da

ADT ast-experimental
Last change on this file since 814a4da was a6bb5fc, checked in by caparsons <caparson@…>, 3 years ago

fixed merge conflict
  • Property mode set to 100644
File size: 11.9 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 // #define SPLITMIX_32
35#else // 32-bit architecture
36 // 64-bit generators
37 //#define XORSHIFT_13_7_17
38 #define XOSHIRO256PP
39
40 // 32-bit generators
41 //#define XORSHIFT_6_21_7
42 #define XOSHIRO128PP
43 // #define SPLITMIX_32
44#endif // __x86_64__
45
46// Define C/CFA PRNG name and random-state.
47
48#ifdef XOSHIRO256PP
49#define PRNG_NAME_64 xoshiro256pp
50#define PRNG_STATE_64_T GLUE(PRNG_NAME_64,_t)
51typedef struct { uint64_t s0, s1, s2, s3; } PRNG_STATE_64_T;
52#endif // XOSHIRO256PP
53
54#ifdef XOSHIRO128PP
55#define PRNG_NAME_32 xoshiro128pp
56#define PRNG_STATE_32_T GLUE(PRNG_NAME_32,_t)
57typedef struct { uint32_t s0, s1, s2, s3; } PRNG_STATE_32_T;
58#endif // XOSHIRO128PP
59
60#ifdef LEHMER64
61#define PRNG_NAME_64 lehmer64
62#define PRNG_STATE_64_T __uint128_t
63#endif // LEHMER64
64
65#ifdef WYHASH64
66#define PRNG_NAME_64 wyhash64
67#define PRNG_STATE_64_T uint64_t
68#endif // LEHMER64
69
70#ifdef XORSHIFT_13_7_17
71#define PRNG_NAME_64 xorshift_13_7_17
72#define PRNG_STATE_64_T uint64_t
73#endif // XORSHIFT_13_7_17
74
75#ifdef XORSHIFT_6_21_7
76#define PRNG_NAME_32 xorshift_6_21_7
77#define PRNG_STATE_32_T uint32_t
78#endif // XORSHIFT_6_21_7
79
80#ifdef SPLITMIX_32
81#define PRNG_NAME_32 splitmix
82#define PRNG_STATE_32_T uint32_t
83#endif // SPLITMIX32
84
85#ifdef XORSHIFT_12_25_27
86#define PRNG_NAME_64 xorshift_12_25_27
87#define PRNG_STATE_64_T uint64_t
88#endif // XORSHIFT_12_25_27
89
90#ifdef KISS_64
91#define PRNG_NAME_64 kiss_64
92#define PRNG_STATE_64_T GLUE(PRNG_NAME_64,_t)
93typedef struct { uint64_t z, w, jsr, jcong; } PRNG_STATE_64_T;
94#endif // KISS_^64
95
96#ifdef XORWOW
97#define PRNG_NAME_32 xorwow
98#define PRNG_STATE_32_T GLUE(PRNG_NAME_32,_t)
99typedef struct { uint32_t a, b, c, d, counter; } PRNG_STATE_32_T;
100#endif // XOSHIRO128PP
101
102#define PRNG_SET_SEED_64 GLUE(PRNG_NAME_64,_set_seed)
103#define PRNG_SET_SEED_32 GLUE(PRNG_NAME_32,_set_seed)
104
105
106// Default PRNG used by runtime.
107#ifdef __x86_64__ // 64-bit architecture
108#define PRNG_NAME PRNG_NAME_64
109#define PRNG_STATE_T PRNG_STATE_64_T
110#else // 32-bit architecture
111#define PRNG_NAME PRNG_NAME_32
112#define PRNG_STATE_T PRNG_STATE_32_T
113#endif // __x86_64__
114
115#define PRNG_SET_SEED GLUE(PRNG_NAME,_set_seed)
116
117
118// ALL PRNG ALGORITHMS ARE OPTIMIZED SO THAT THE PRNG LOGIC CAN HAPPEN IN PARALLEL WITH THE USE OF THE RESULT.
119// Specifically, the current random state is copied for returning, before computing the next value. As a consequence,
120// the set_seed routine primes the PRNG by calling it with the state so the seed is not return as the first random
121// value.
122
123
124#ifdef __cforall // don't include in C code (invoke.h)
125
126// Splitmix32
127// https://github.com/bryc/code/blob/master/jshash/PRNGs.md#splitmix32
128// Splitmix32 is not recommended for demanding random number requirements,
129// but is often used to calculate initial states for other more complex
130// pseudo-random number generators.
131// SplitMix32 is a 32 bit variant of Splitmix64
132
133static inline uint32_t splitmix32( uint32_t & state ) {
134 state += 0x9e3779b9;
135 uint64_t z = state;
136 z = (z ^ (z >> 15)) * 0x85ebca6b;
137 z = (z ^ (z >> 13)) * 0xc2b2ae35;
138 return z ^ (z >> 16);
139}
140
141static inline void splitmix32_set_seed( uint32_t & state , uint64_t seed ) {
142 state = seed;
143 splitmix32( state ); // prime
144} // splitmix32_set_seed
145
146#ifdef __SIZEOF_INT128__
147 //--------------------------------------------------
148 static inline uint64_t lehmer64( __uint128_t & state ) {
149 __uint128_t ret = state;
150 state *= 0x_da94_2042_e4dd_58b5;
151 return ret >> 64;
152 } // lehmer64
153
154 static inline void lehmer64_set_seed( __uint128_t & state, uint64_t seed ) {
155 // The seed needs to be coprime with the 2^64 modulus to get the largest period, so no factors of 2 in the seed.
156 state = seed;
157 lehmer64( state ); // prime
158 } // lehmer64_set_seed
159
160 //--------------------------------------------------
161 static inline uint64_t wyhash64( uint64_t & state ) {
162 uint64_t ret = state;
163 state += 0x_60be_e2be_e120_fc15;
164 __uint128_t tmp;
165 tmp = (__uint128_t) ret * 0x_a3b1_9535_4a39_b70d;
166 uint64_t m1 = (tmp >> 64) ^ tmp;
167 tmp = (__uint128_t)m1 * 0x_1b03_7387_12fa_d5c9;
168 uint64_t m2 = (tmp >> 64) ^ tmp;
169 return m2;
170 } // wyhash64
171
172 static inline void wyhash64_set_seed( uint64_t & state, uint64_t seed ) {
173 state = seed;
174 wyhash64( state ); // prime
175 } // wyhash64_set_seed
176#endif // __SIZEOF_INT128__
177
178// https://prng.di.unimi.it/xoshiro256starstar.c
179//
180// This is xoshiro256++ 1.0, one of our all-purpose, rock-solid generators. It has excellent (sub-ns) speed, a state
181// (256 bits) that is large enough for any parallel application, and it passes all tests we are aware of.
182//
183// For generating just floating-point numbers, xoshiro256+ is even faster.
184//
185// The state must be seeded so that it is not everywhere zero. If you have a 64-bit seed, we suggest to seed a
186// splitmix64 generator and use its output to fill s.
187
188#ifndef XOSHIRO256PP
189typedef struct { uint64_t s0, s1, s2, s3; } xoshiro256pp_t;
190#endif // ! XOSHIRO256PP
191
192static inline uint64_t xoshiro256pp( xoshiro256pp_t & rs ) with(rs) {
193 inline uint64_t rotl( const uint64_t x, int k ) {
194 return (x << k) | (x >> (64 - k));
195 } // rotl
196
197 const uint64_t result = rotl( s0 + s3, 23 ) + s0;
198 const uint64_t t = s1 << 17;
199
200 s2 ^= s0;
201 s3 ^= s1;
202 s1 ^= s2;
203 s0 ^= s3;
204 s2 ^= t;
205 s3 = rotl( s3, 45 );
206 return result;
207} // xoshiro256pp
208
209static inline void xoshiro256pp_set_seed( xoshiro256pp_t & state, uint64_t seed ) {
210 uint64_t state;
211 wyhash64_set_seed( state, seed );
212 // these are done explicitly in this order to attain repeatable seeding.
213 // do not call splitmix32 directly in the state init since order of argument evaluation
214 // may not be consistent leading to irreproducible seeding
215 uint64_t seed1 = wyhash64( state );
216 uint64_t seed2 = wyhash64( state );
217 uint64_t seed3 = wyhash64( state );
218 uint64_t seed4 = wyhash64( state );
219 state = (xoshiro256pp_t){ seed1, seed2, seed3, seed4 };
220 xoshiro256pp( state ); // prime
221} // xoshiro256pp_set_seed
222
223// https://prng.di.unimi.it/xoshiro128plusplus.c
224//
225// This is xoshiro128++ 1.0, one of our 32-bit all-purpose, rock-solid generators. It has excellent speed, a state size
226// (128 bits) that is large enough for mild parallelism, and it passes all tests we are aware of.
227//
228// For generating just single-precision (i.e., 32-bit) floating-point numbers, xoshiro128+ is even faster.
229//
230// The state must be seeded so that it is not everywhere zero.
231
232#ifndef XOSHIRO128PP
233typedef struct { uint32_t s0, s1, s2, s3; } xoshiro128pp_t;
234#endif // ! XOSHIRO128PP
235
236static inline uint32_t xoshiro128pp( xoshiro128pp_t & rs ) with(rs) {
237 inline uint32_t rotl( const uint32_t x, int k ) {
238 return (x << k) | (x >> (32 - k));
239 } // rotl
240
241 const uint32_t result = rotl( s0 + s3, 7 ) + s0;
242 const uint32_t t = s1 << 9;
243
244 s2 ^= s0;
245 s3 ^= s1;
246 s1 ^= s2;
247 s0 ^= s3;
248 s2 ^= t;
249 s3 = rotl( s3, 11 );
250 return result;
251} // xoshiro128pp
252
253static inline void xoshiro128pp_set_seed( xoshiro128pp_t & state, uint32_t seed ) {
254 // these are done explicitly in this order to attain repeatable seeding.
255 // do not call splitmix32 directly in the state init since order of argument evaluation
256 // may not be consistent leading to irreproducible seeding
257 uint32_t seed1 = splitmix32( seed );
258 uint32_t seed2 = splitmix32( seed );
259 uint32_t seed3 = splitmix32( seed );
260 uint32_t seed4 = splitmix32( seed );
261 state = (xoshiro128pp_t){ seed1, seed2, seed3, seed4 };
262 xoshiro128pp( state ); // prime
263} // xoshiro128pp_set_seed
264
265//--------------------------------------------------
266static inline uint64_t xorshift_13_7_17( uint64_t & state ) {
267 uint64_t ret = state;
268 state ^= state << 13;
269 state ^= state >> 7;
270 state ^= state << 17;
271 return ret;
272} // xorshift_13_7_17
273
274static inline void xorshift_13_7_17_set_seed( uint64_t & state, uint64_t seed ) {
275 state = seed;
276 xorshift_13_7_17( state ); // prime
277} // xorshift_13_7_17_set_seed
278
279//--------------------------------------------------
280// Marsaglia shift-XOR PRNG with thread-local state
281// Period is 4G-1
282// 0 is absorbing and must be avoided
283// Low-order bits are not particularly random
284static inline uint32_t xorshift_6_21_7( uint32_t & state ) {
285 uint32_t ret = state;
286 state ^= state << 6;
287 state ^= state >> 21;
288 state ^= state << 7;
289 return ret;
290} // xorshift_6_21_7
291
292static inline void xorshift_6_21_7_set_seed( uint32_t & state, uint32_t seed ) {
293 state = seed;
294 xorshift_6_21_7( state ); // prime
295} // xorshift_6_21_7_set_seed
296
297//--------------------------------------------------
298// The state must be seeded with a nonzero value.
299static inline uint64_t xorshift_12_25_27( uint64_t & state ) {
300 uint64_t ret = state;
301 state ^= state >> 12;
302 state ^= state << 25;
303 state ^= state >> 27;
304 return ret * 0x_2545_F491_4F6C_DD1D;
305} // xorshift_12_25_27
306
307static inline void xorshift_12_25_27_set_seed( uint64_t & state, uint64_t seed ) {
308 state = seed;
309 xorshift_12_25_27( state ); // prime
310} // xorshift_12_25_27_set_seed
311
312//--------------------------------------------------
313// The state must be seeded with a nonzero value.
314#ifndef KISS_64
315typedef struct { uint64_t z, w, jsr, jcong; } kiss_64_t;
316#endif // ! KISS_64
317
318static inline uint64_t kiss_64( kiss_64_t & rs ) with(rs) {
319 kiss_64_t ret = rs;
320 z = 36969 * (z & 65535) + (z >> 16);
321 w = 18000 * (w & 65535) + (w >> 16);
322 jsr ^= (jsr << 13);
323 jsr ^= (jsr >> 17);
324 jsr ^= (jsr << 5);
325 jcong = 69069 * jcong + 1234567;
326 return (((ret.z << 16) + ret.w) ^ ret.jcong) + ret.jsr;
327} // kiss_64
328
329static inline void kiss_64_set_seed( kiss_64_t & rs, uint64_t seed ) with(rs) {
330 z = 1; w = 1; jsr = 4; jcong = seed;
331 kiss_64( rs ); // prime
332} // kiss_64_set_seed
333
334//--------------------------------------------------
335// The state array must be initialized to non-zero in the first four words.
336#ifndef XORWOW
337typedef struct { uint32_t a, b, c, d, counter; } xorwow_t;
338#endif // ! XORWOW
339
340static inline uint32_t xorwow( xorwow_t & rs ) with(rs) {
341 // Algorithm "xorwow" from p. 5 of Marsaglia, "Xorshift RNGs".
342 uint32_t ret = a + counter;
343 uint32_t t = d;
344
345 uint32_t const s = a;
346 d = c;
347 c = b;
348 b = s;
349
350 t ^= t >> 2;
351 t ^= t << 1;
352 t ^= s ^ (s << 4);
353 a = t;
354 counter += 362437;
355 return ret;
356} // xorwow
357
358static inline void xorwow_set_seed( xorwow_t & rs, uint32_t seed ) {
359 // these are done explicitly in this order to attain repeatable seeding.
360 // do not call splitmix32 directly in the state init since order of argument evaluation
361 // may not be consistent leading to irreproducible seeding
362 uint32_t seed1 = splitmix32( seed );
363 uint32_t seed2 = splitmix32( seed );
364 uint32_t seed3 = splitmix32( seed );
365 uint32_t seed4 = splitmix32( seed );
366 rs = (xorwow_t){ seed1, seed2, seed3, seed4, 0 };
367 xorwow( rs ); // prime
368} // xorwow_set_seed
369
370//--------------------------------------------------
371// Used in __tls_rand_fwd
372#define M (1_l64u << 48_l64u)
373#define A (25214903917_l64u)
374#define AI (18446708753438544741_l64u)
375#define C (11_l64u)
376#define D (16_l64u)
377
378// Bi-directional LCG random-number generator
379static inline uint32_t LCGBI_fwd( uint64_t & rs ) {
380 rs = (A * rs + C) & (M - 1);
381 return rs >> D;
382} // LCGBI_fwd
383
384static inline uint32_t LCGBI_bck( uint64_t & rs ) {
385 unsigned int r = rs >> D;
386 rs = AI * (rs - C) & (M - 1);
387 return r;
388} // LCGBI_bck
389
390#undef M
391#undef A
392#undef AI
393#undef C
394#undef D
395
396#endif // __cforall
Note: See TracBrowser for help on using the repository browser.