source: libcfa/src/bits/random.hfa @ 4c6ba5a

ADTast-experimental
Last change on this file since 4c6ba5a was 4c6ba5a, checked in by caparsons <caparson@…>, 8 months ago

refactored to use generators for seeding state, added splitmix32 for 32 bit seeding

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