source: libcfa/src/bits/random.hfa @ 90fb672

ADTast-experimental
Last change on this file since 90fb672 was 90fb672, checked in by Peter A. Buhr <pabuhr@…>, 19 months ago

use splitmix32/64 to prime set seed for all PRNG

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