1
|
|
/*
|
2
|
|
* PCG64 Random Number Generation for C.
|
3
|
|
*
|
4
|
|
* Copyright 2014 Melissa O'Neill <oneill@pcg-random.org>
|
5
|
|
* Copyright 2015 Robert Kern <robert.kern@gmail.com>
|
6
|
|
*
|
7
|
|
* Licensed under the Apache License, Version 2.0 (the "License");
|
8
|
|
* you may not use this file except in compliance with the License.
|
9
|
|
* You may obtain a copy of the License at
|
10
|
|
*
|
11
|
|
* http://www.apache.org/licenses/LICENSE-2.0
|
12
|
|
*
|
13
|
|
* Unless required by applicable law or agreed to in writing, software
|
14
|
|
* distributed under the License is distributed on an "AS IS" BASIS,
|
15
|
|
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
16
|
|
* See the License for the specific language governing permissions and
|
17
|
|
* limitations under the License.
|
18
|
|
*
|
19
|
|
* For additional information about the PCG random number generation scheme,
|
20
|
|
* including its license and other licensing options, visit
|
21
|
|
*
|
22
|
|
* http://www.pcg-random.org
|
23
|
|
*
|
24
|
|
* Relicensed MIT in May 2019
|
25
|
|
*
|
26
|
|
* The MIT License
|
27
|
|
*
|
28
|
|
* PCG Random Number Generation for C.
|
29
|
|
*
|
30
|
|
* Copyright 2014 Melissa O'Neill <oneill@pcg-random.org>
|
31
|
|
*
|
32
|
|
* Permission is hereby granted, free of charge, to any person obtaining
|
33
|
|
* a copy of this software and associated documentation files (the "Software"),
|
34
|
|
* to deal in the Software without restriction, including without limitation
|
35
|
|
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
|
36
|
|
* and/or sell copies of the Software, and to permit persons to whom the
|
37
|
|
* Software is furnished to do so, subject to the following conditions:
|
38
|
|
*
|
39
|
|
* The above copyright notice and this permission notice shall be included in
|
40
|
|
* all copies or substantial portions of the Software.
|
41
|
|
*
|
42
|
|
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
43
|
|
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
|
44
|
|
* FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
|
45
|
|
* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
|
46
|
|
* IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
|
47
|
|
* CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
|
48
|
|
*/
|
49
|
|
|
50
|
|
#ifndef PCG64_H_INCLUDED
|
51
|
|
#define PCG64_H_INCLUDED 1
|
52
|
|
|
53
|
|
#include <inttypes.h>
|
54
|
|
|
55
|
|
#ifdef _WIN32
|
56
|
|
#include <stdlib.h>
|
57
|
|
#define inline __forceinline
|
58
|
|
#endif
|
59
|
|
|
60
|
|
#if __GNUC_GNU_INLINE__ && !defined(__cplusplus)
|
61
|
|
#error Nonstandard GNU inlining semantics. Compile with -std=c99 or better.
|
62
|
|
#endif
|
63
|
|
|
64
|
|
#if __cplusplus
|
65
|
|
extern "C" {
|
66
|
|
#endif
|
67
|
|
|
68
|
|
#if defined(__SIZEOF_INT128__) && !defined(PCG_FORCE_EMULATED_128BIT_MATH)
|
69
|
|
typedef __uint128_t pcg128_t;
|
70
|
|
#define PCG_128BIT_CONSTANT(high, low) (((pcg128_t)(high) << 64) + low)
|
71
|
|
#else
|
72
|
|
typedef struct {
|
73
|
|
uint64_t high;
|
74
|
|
uint64_t low;
|
75
|
|
} pcg128_t;
|
76
|
|
|
77
|
|
static inline pcg128_t PCG_128BIT_CONSTANT(uint64_t high, uint64_t low) {
|
78
|
|
pcg128_t result;
|
79
|
|
result.high = high;
|
80
|
|
result.low = low;
|
81
|
|
return result;
|
82
|
|
}
|
83
|
|
|
84
|
|
#define PCG_EMULATED_128BIT_MATH 1
|
85
|
|
#endif
|
86
|
|
|
87
|
|
typedef struct { pcg128_t state; } pcg_state_128;
|
88
|
|
|
89
|
|
typedef struct {
|
90
|
|
pcg128_t state;
|
91
|
|
pcg128_t inc;
|
92
|
|
} pcg_state_setseq_128;
|
93
|
|
|
94
|
|
#define PCG_DEFAULT_MULTIPLIER_HIGH 2549297995355413924ULL
|
95
|
|
#define PCG_DEFAULT_MULTIPLIER_LOW 4865540595714422341ULL
|
96
|
|
|
97
|
|
#define PCG_DEFAULT_MULTIPLIER_128 \
|
98
|
|
PCG_128BIT_CONSTANT(PCG_DEFAULT_MULTIPLIER_HIGH, PCG_DEFAULT_MULTIPLIER_LOW)
|
99
|
|
#define PCG_DEFAULT_INCREMENT_128 \
|
100
|
|
PCG_128BIT_CONSTANT(6364136223846793005ULL, 1442695040888963407ULL)
|
101
|
|
#define PCG_STATE_SETSEQ_128_INITIALIZER \
|
102
|
|
{ \
|
103
|
|
PCG_128BIT_CONSTANT(0x979c9a98d8462005ULL, 0x7d3e9cb6cfe0549bULL) \
|
104
|
|
, PCG_128BIT_CONSTANT(0x0000000000000001ULL, 0xda3e39cb94b95bdbULL) \
|
105
|
|
}
|
106
|
|
|
107
|
|
static inline uint64_t pcg_rotr_64(uint64_t value, unsigned int rot) {
|
108
|
|
#ifdef _WIN32
|
109
|
|
return _rotr64(value, rot);
|
110
|
|
#else
|
111
|
1
|
return (value >> rot) | (value << ((-rot) & 63));
|
112
|
|
#endif
|
113
|
|
}
|
114
|
|
|
115
|
|
#ifdef PCG_EMULATED_128BIT_MATH
|
116
|
|
|
117
|
|
static inline pcg128_t pcg128_add(pcg128_t a, pcg128_t b) {
|
118
|
|
pcg128_t result;
|
119
|
|
|
120
|
|
result.low = a.low + b.low;
|
121
|
|
result.high = a.high + b.high + (result.low < b.low);
|
122
|
|
return result;
|
123
|
|
}
|
124
|
|
|
125
|
|
static inline void _pcg_mult64(uint64_t x, uint64_t y, uint64_t *z1,
|
126
|
|
uint64_t *z0) {
|
127
|
|
|
128
|
|
#if defined _WIN32 && _MSC_VER >= 1900 && _M_AMD64
|
129
|
|
z0[0] = _umul128(x, y, z1);
|
130
|
|
#else
|
131
|
|
uint64_t x0, x1, y0, y1;
|
132
|
|
uint64_t w0, w1, w2, t;
|
133
|
|
/* Lower 64 bits are straightforward clock-arithmetic. */
|
134
|
|
*z0 = x * y;
|
135
|
|
|
136
|
|
x0 = x & 0xFFFFFFFFULL;
|
137
|
|
x1 = x >> 32;
|
138
|
|
y0 = y & 0xFFFFFFFFULL;
|
139
|
|
y1 = y >> 32;
|
140
|
|
w0 = x0 * y0;
|
141
|
|
t = x1 * y0 + (w0 >> 32);
|
142
|
|
w1 = t & 0xFFFFFFFFULL;
|
143
|
|
w2 = t >> 32;
|
144
|
|
w1 += x0 * y1;
|
145
|
|
*z1 = x1 * y1 + w2 + (w1 >> 32);
|
146
|
|
#endif
|
147
|
|
}
|
148
|
|
|
149
|
|
static inline pcg128_t pcg128_mult(pcg128_t a, pcg128_t b) {
|
150
|
|
uint64_t h1;
|
151
|
|
pcg128_t result;
|
152
|
|
|
153
|
|
h1 = a.high * b.low + a.low * b.high;
|
154
|
|
_pcg_mult64(a.low, b.low, &(result.high), &(result.low));
|
155
|
|
result.high += h1;
|
156
|
|
return result;
|
157
|
|
}
|
158
|
|
|
159
|
|
static inline void pcg_setseq_128_step_r(pcg_state_setseq_128 *rng) {
|
160
|
|
rng->state = pcg128_add(pcg128_mult(rng->state, PCG_DEFAULT_MULTIPLIER_128),
|
161
|
|
rng->inc);
|
162
|
|
}
|
163
|
|
|
164
|
|
static inline uint64_t pcg_output_xsl_rr_128_64(pcg128_t state) {
|
165
|
|
return pcg_rotr_64(state.high ^ state.low, state.high >> 58u);
|
166
|
|
}
|
167
|
|
|
168
|
|
static inline void pcg_setseq_128_srandom_r(pcg_state_setseq_128 *rng,
|
169
|
|
pcg128_t initstate,
|
170
|
|
pcg128_t initseq) {
|
171
|
|
rng->state = PCG_128BIT_CONSTANT(0ULL, 0ULL);
|
172
|
|
rng->inc.high = initseq.high << 1u;
|
173
|
|
rng->inc.high |= initseq.low >> 63u;
|
174
|
|
rng->inc.low = (initseq.low << 1u) | 1u;
|
175
|
|
pcg_setseq_128_step_r(rng);
|
176
|
|
rng->state = pcg128_add(rng->state, initstate);
|
177
|
|
pcg_setseq_128_step_r(rng);
|
178
|
|
}
|
179
|
|
|
180
|
|
static inline uint64_t
|
181
|
|
pcg_setseq_128_xsl_rr_64_random_r(pcg_state_setseq_128 *rng) {
|
182
|
|
#if defined _WIN32 && _MSC_VER >= 1900 && _M_AMD64
|
183
|
|
uint64_t h1;
|
184
|
|
pcg128_t product;
|
185
|
|
|
186
|
|
/* Manually inline the multiplication and addition using intrinsics */
|
187
|
|
h1 = rng->state.high * PCG_DEFAULT_MULTIPLIER_LOW +
|
188
|
|
rng->state.low * PCG_DEFAULT_MULTIPLIER_HIGH;
|
189
|
|
product.low =
|
190
|
|
_umul128(rng->state.low, PCG_DEFAULT_MULTIPLIER_LOW, &(product.high));
|
191
|
|
product.high += h1;
|
192
|
|
_addcarry_u64(_addcarry_u64(0, product.low, rng->inc.low, &(rng->state.low)),
|
193
|
|
product.high, rng->inc.high, &(rng->state.high));
|
194
|
|
return _rotr64(rng->state.high ^ rng->state.low, rng->state.high >> 58u);
|
195
|
|
#else
|
196
|
|
pcg_setseq_128_step_r(rng);
|
197
|
|
return pcg_output_xsl_rr_128_64(rng->state);
|
198
|
|
#endif
|
199
|
|
}
|
200
|
|
|
201
|
|
#else /* PCG_EMULATED_128BIT_MATH */
|
202
|
|
|
203
|
|
static inline void pcg_setseq_128_step_r(pcg_state_setseq_128 *rng) {
|
204
|
1
|
rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_128 + rng->inc;
|
205
|
|
}
|
206
|
|
|
207
|
|
static inline uint64_t pcg_output_xsl_rr_128_64(pcg128_t state) {
|
208
|
1
|
return pcg_rotr_64(((uint64_t)(state >> 64u)) ^ (uint64_t)state,
|
209
|
1
|
state >> 122u);
|
210
|
|
}
|
211
|
|
|
212
|
|
static inline uint64_t
|
213
|
|
pcg_setseq_128_xsl_rr_64_random_r(pcg_state_setseq_128* rng)
|
214
|
|
{
|
215
|
1
|
pcg_setseq_128_step_r(rng);
|
216
|
1
|
return pcg_output_xsl_rr_128_64(rng->state);
|
217
|
|
}
|
218
|
|
|
219
|
|
static inline void pcg_setseq_128_srandom_r(pcg_state_setseq_128 *rng,
|
220
|
|
pcg128_t initstate,
|
221
|
|
pcg128_t initseq) {
|
222
|
|
rng->state = 0U;
|
223
|
1
|
rng->inc = (initseq << 1u) | 1u;
|
224
|
1
|
pcg_setseq_128_step_r(rng);
|
225
|
1
|
rng->state += initstate;
|
226
|
1
|
pcg_setseq_128_step_r(rng);
|
227
|
|
}
|
228
|
|
|
229
|
|
#endif /* PCG_EMULATED_128BIT_MATH */
|
230
|
|
|
231
|
|
static inline uint64_t
|
232
|
|
pcg_setseq_128_xsl_rr_64_boundedrand_r(pcg_state_setseq_128 *rng,
|
233
|
|
uint64_t bound) {
|
234
|
|
uint64_t threshold = -bound % bound;
|
235
|
|
for (;;) {
|
236
|
|
uint64_t r = pcg_setseq_128_xsl_rr_64_random_r(rng);
|
237
|
|
if (r >= threshold)
|
238
|
|
return r % bound;
|
239
|
|
}
|
240
|
|
}
|
241
|
|
|
242
|
|
extern pcg128_t pcg_advance_lcg_128(pcg128_t state, pcg128_t delta,
|
243
|
|
pcg128_t cur_mult, pcg128_t cur_plus);
|
244
|
|
|
245
|
|
static inline void pcg_setseq_128_advance_r(pcg_state_setseq_128 *rng,
|
246
|
|
pcg128_t delta) {
|
247
|
1
|
rng->state = pcg_advance_lcg_128(rng->state, delta,
|
248
|
|
PCG_DEFAULT_MULTIPLIER_128, rng->inc);
|
249
|
|
}
|
250
|
|
|
251
|
|
typedef pcg_state_setseq_128 pcg64_random_t;
|
252
|
|
#define pcg64_random_r pcg_setseq_128_xsl_rr_64_random_r
|
253
|
|
#define pcg64_boundedrand_r pcg_setseq_128_xsl_rr_64_boundedrand_r
|
254
|
|
#define pcg64_srandom_r pcg_setseq_128_srandom_r
|
255
|
|
#define pcg64_advance_r pcg_setseq_128_advance_r
|
256
|
|
#define PCG64_INITIALIZER PCG_STATE_SETSEQ_128_INITIALIZER
|
257
|
|
|
258
|
|
#if __cplusplus
|
259
|
|
}
|
260
|
|
#endif
|
261
|
|
|
262
|
|
typedef struct s_pcg64_state {
|
263
|
|
pcg64_random_t *pcg_state;
|
264
|
|
int has_uint32;
|
265
|
|
uint32_t uinteger;
|
266
|
|
} pcg64_state;
|
267
|
|
|
268
|
|
static inline uint64_t pcg64_next64(pcg64_state *state) {
|
269
|
1
|
return pcg64_random_r(state->pcg_state);
|
270
|
|
}
|
271
|
|
|
272
|
|
static inline uint32_t pcg64_next32(pcg64_state *state) {
|
273
|
|
uint64_t next;
|
274
|
1
|
if (state->has_uint32) {
|
275
|
1
|
state->has_uint32 = 0;
|
276
|
1
|
return state->uinteger;
|
277
|
|
}
|
278
|
1
|
next = pcg64_random_r(state->pcg_state);
|
279
|
1
|
state->has_uint32 = 1;
|
280
|
1
|
state->uinteger = (uint32_t)(next >> 32);
|
281
|
1
|
return (uint32_t)(next & 0xffffffff);
|
282
|
|
}
|
283
|
|
|
284
|
|
void pcg64_advance(pcg64_state *state, uint64_t *step);
|
285
|
|
|
286
|
|
void pcg64_set_seed(pcg64_state *state, uint64_t *seed, uint64_t *inc);
|
287
|
|
|
288
|
|
void pcg64_get_state(pcg64_state *state, uint64_t *state_arr, int *has_uint32,
|
289
|
|
uint32_t *uinteger);
|
290
|
|
|
291
|
|
void pcg64_set_state(pcg64_state *state, uint64_t *state_arr, int has_uint32,
|
292
|
|
uint32_t uinteger);
|
293
|
|
|
294
|
|
#endif /* PCG64_H_INCLUDED */
|