Context Navigation

random.hfa @ cd477ca

ADTast-experimental

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format