numpy / numpy
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 */

Read our documentation on viewing source code .

Loading