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

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

remove 32-bit use of lehmer64, which requires uint128_t

  • Property mode set to 100644
File size: 6.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 : Thu Dec  1 11:52:08 2022
13// Update Count     : 115
14//
15
16#pragma once
17
18#include <stdint.h>
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#define LEHMER64
26#define XORSHIFT_6_21_7
27//#define XOSHIRO256PP
28//#define XOSHIRO128PP
29#else                                                                                                   // 32-bit architecture
30#define XORSHIFT_13_7_17
31#define XORSHIFT_6_21_7
32//#define XOSHIRO256PP
33//#define XOSHIRO128PP
34#endif // __x86_64__
35
36// C/CFA PRNG name and random-state.
37
38#ifdef LEHMER64
39#define PRNG_NAME_64 lehmer64
40#define PRNG_STATE_64_T __uint128_t
41#endif // LEHMER64
42
43#ifdef XORSHIFT_13_7_17
44#define PRNG_NAME_64 xorshift_13_7_17
45#define PRNG_STATE_64_T uint64_t
46#endif // XORSHIFT_13_7_17
47
48#ifdef XORSHIFT_6_21_7
49#define PRNG_NAME_32 xorshift_6_21_7
50#define PRNG_STATE_32_T uint32_t
51#endif // XORSHIFT_6_21_7
52
53#ifdef XOSHIRO256PP
54#define PRNG_NAME_64 xoshiro256pp
55#define PRNG_STATE_64_T struct GLUE(PRNG_NAME_64,_t)
56PRNG_STATE_64_T { uint64_t s[4]; };
57#endif // XOSHIRO256PP
58
59#ifdef XOSHIRO128PP
60#define PRNG_NAME_32 xoshiro128pp
61#define PRNG_STATE_32_T struct GLUE(PRNG_NAME_32,_t)
62PRNG_STATE_32_T { uint32_t s[4]; };
63#endif // XOSHIRO128PP
64
65#define PRNG_SET_SEED_64 GLUE(PRNG_NAME_64,_set_seed)
66#define PRNG_SET_SEED_32 GLUE(PRNG_NAME_32,_set_seed)
67
68
69// Default PRNG used by runtime.
70#ifdef __x86_64__                                                                               // 64-bit architecture
71#define PRNG_NAME PRNG_NAME_64
72#define PRNG_STATE_T PRNG_STATE_64_T
73#else                                                                                                   // 32-bit architecture
74#define PRNG_NAME PRNG_NAME_32
75#define PRNG_STATE_T PRNG_STATE_32_T
76#endif // __x86_64__
77
78#define PRNG_SET_SEED GLUE(PRNG_NAME,_set_seed)
79
80
81#ifdef __cforall                                                                                // don't include in C code (invoke.h)
82
83// https://prng.di.unimi.it/xoshiro128plusplus.c
84//
85// This is xoshiro128++ 1.0, one of our 32-bit all-purpose, rock-solid generators. It has excellent speed, a state size
86// (128 bits) that is large enough for mild parallelism, and it passes all tests we are aware of.
87//
88// For generating just single-precision (i.e., 32-bit) floating-point numbers, xoshiro128+ is even faster.
89//
90// The state must be seeded so that it is not everywhere zero.
91
92#ifndef XOSHIRO128PP
93struct xoshiro128pp_t { uint32_t s[4]; };
94#endif // ! XOSHIRO128PP
95
96static inline uint32_t xoshiro128pp( xoshiro128pp_t & rs ) with(rs) {
97        inline uint32_t rotl( const uint32_t x, int k ) {
98                return (x << k) | (x >> (32 - k));
99        }
100
101        const uint32_t result = rotl( s[0] + s[3], 7 ) + s[0];
102        const uint32_t t = s[1] << 9;
103
104        s[2] ^= s[0];
105        s[3] ^= s[1];
106        s[1] ^= s[2];
107        s[0] ^= s[3];
108        s[2] ^= t;
109        s[3] = rotl( s[3], 11 );
110        return result;
111}
112
113static inline void xoshiro128pp_set_seed( xoshiro128pp_t & state, uint32_t seed ) {
114        state = (xoshiro128pp_t){ {seed, seed, seed, seed} };
115} // xoshiro128pp_set_seed
116
117// This is xoshiro256++ 1.0, one of our all-purpose, rock-solid generators.  It has excellent (sub-ns) speed, a state
118// (256 bits) that is large enough for any parallel application, and it passes all tests we are aware of.
119//
120// For generating just floating-point numbers, xoshiro256+ is even faster.
121//
122// The state must be seeded so that it is not everywhere zero. If you have a 64-bit seed, we suggest to seed a
123// splitmix64 generator and use its output to fill s.
124
125#ifndef XOSHIRO256PP
126struct xoshiro256pp_t { uint64_t s[4]; };
127#endif // ! XOSHIRO256PP
128
129static inline uint64_t xoshiro256pp( xoshiro256pp_t & rs ) with(rs) {
130        inline uint64_t rotl(const uint64_t x, int k) {
131                return (x << k) | (x >> (64 - k));
132        }
133
134        const uint64_t result = rotl( s[0] + s[3], 23 ) + s[0];
135        const uint64_t t = s[1] << 17;
136
137        s[2] ^= s[0];
138        s[3] ^= s[1];
139        s[1] ^= s[2];
140        s[0] ^= s[3];
141        s[2] ^= t;
142        s[3] = rotl( s[3], 45 );
143        return result;
144}
145
146static inline void xoshiro256pp_set_seed( xoshiro256pp_t & state,  uint64_t seed ) {
147        state = (xoshiro256pp_t){ {seed, seed, seed, seed} };
148} // xoshiro256pp_set_seed
149
150#ifdef __SIZEOF_INT128__
151        // Pipelined to allow out-of-order overlap with reduced dependencies. Critically, the current random state is
152        // returned (copied), and then compute and store the next random value.
153        //--------------------------------------------------
154        static inline uint64_t lehmer64( __uint128_t & state ) {
155                __uint128_t ret = state;
156                state *= 0xda942042e4dd58b5;
157                return ret >> 64;
158        } // lehmer64
159
160        static inline void lehmer64_set_seed( __uint128_t & state, uint64_t seed ) {
161                state = seed;
162        } // lehmer64_set_seed
163
164        //--------------------------------------------------
165        static inline uint64_t wyhash64( uint64_t & state ) {
166                state += 0x60bee2bee120fc15;
167                __uint128_t tmp;
168                tmp = (__uint128_t) state * 0xa3b195354a39b70d;
169                uint64_t m1 = (tmp >> 64) ^ tmp;
170                tmp = (__uint128_t)m1 * 0x1b03738712fad5c9;
171                uint64_t m2 = (tmp >> 64) ^ tmp;
172                return m2;
173        }
174
175        static inline void wyhash64_set_seed( __uint128_t & state, uint64_t seed ) {
176                state = seed;
177        } // lehmer64_set_seed
178#endif // __SIZEOF_INT128__
179
180//--------------------------------------------------
181static inline uint64_t xorshift_13_7_17( uint64_t & state ) {
182        uint64_t ret = state;
183        state ^= state << 13;
184        state ^= state >> 7;
185        state ^= state << 17;
186        return ret;
187}
188
189static inline void xorshift_13_7_17_set_seed( uint64_t & state, uint32_t seed ) {
190        state = seed;
191}
192
193//--------------------------------------------------
194static inline uint32_t xorshift_6_21_7( uint32_t & state ) {
195        uint32_t ret = state;
196        state ^= state << 6;
197        state ^= state >> 21;
198        state ^= state << 7;
199        return ret;
200} // xorshift_6_21_7
201
202static inline void xorshift_6_21_7_set_seed( uint32_t & state, uint32_t seed ) {
203        state = seed;
204}
205
206//--------------------------------------------------
207typedef struct {
208        uint32_t a, b, c, d;
209        uint32_t counter;
210} xorwow__state_t;
211
212// The state array must be initialized to not be all zero in the first four words.
213static inline uint32_t xorwow( xorwow__state_t & state ) {
214        // Algorithm "xorwow" from p. 5 of Marsaglia, "Xorshift RNGs".
215        uint32_t ret = state.a + state.counter;
216        uint32_t t = state.d;
217
218        uint32_t const s = state.a;
219        state.d = state.c;
220        state.c = state.b;
221        state.b = s;
222
223        t ^= t >> 2;
224        t ^= t << 1;
225        t ^= s ^ (s << 4);
226        state.a = t;
227
228        state.counter += 362437;
229        return ret;
230}
231
232// Used in __tls_rand_fwd
233//--------------------------------------------------
234#define M  (1_l64u << 48_l64u)
235#define A  (25214903917_l64u)
236#define AI (18446708753438544741_l64u)
237#define C  (11_l64u)
238#define D  (16_l64u)
239
240// Bi-directional LCG random-number generator
241static inline uint32_t LCGBI_fwd( uint64_t & state ) {
242        state = (A * state + C) & (M - 1);
243        return state >> D;
244}
245
246static inline uint32_t LCGBI_bck( uint64_t & state ) {
247        unsigned int r = state >> D;
248        state = AI * (state - C) & (M - 1);
249        return r;
250}
251
252#undef M
253#undef A
254#undef AI
255#undef C
256#undef D
257
258#endif // __cforall
Note: See TracBrowser for help on using the repository browser.