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

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

generalization of PRNG

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