pcg64.h at 5bec8df ⋅ numpy/numpy

/*
 * PCG64 Random Number Generation for C.
 *
 * Copyright 2014 Melissa O'Neill <oneill@pcg-random.org>
 * Copyright 2015 Robert Kern <robert.kern@gmail.com>
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 * For additional information about the PCG random number generation scheme,
 * including its license and other licensing options, visit
 *
 *     http://www.pcg-random.org
 *
 * Relicensed MIT in May 2019
 *
 * The MIT License
 *
 * PCG Random Number Generation for C.
 *
 * Copyright 2014 Melissa O'Neill <oneill@pcg-random.org>
 *
 * Permission is hereby granted, free of charge, to any person obtaining
 * a copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
 * FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
 * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
 * IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
 * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */

#ifndef PCG64_H_INCLUDED
#define PCG64_H_INCLUDED 1

#include <inttypes.h>

#ifdef _WIN32
#include <stdlib.h>
#define inline __forceinline
#endif

#if __GNUC_GNU_INLINE__ && !defined(__cplusplus)
#error Nonstandard GNU inlining semantics. Compile with -std=c99 or better.
#endif

#if __cplusplus
extern "C" {
#endif

#if defined(__SIZEOF_INT128__) && !defined(PCG_FORCE_EMULATED_128BIT_MATH)
typedef __uint128_t pcg128_t;
#define PCG_128BIT_CONSTANT(high, low) (((pcg128_t)(high) << 64) + low)
#else
typedef struct {
  uint64_t high;
  uint64_t low;
} pcg128_t;

static inline pcg128_t PCG_128BIT_CONSTANT(uint64_t high, uint64_t low) {
  pcg128_t result;
  result.high = high;
  result.low = low;
  return result;
}

#define PCG_EMULATED_128BIT_MATH 1
#endif

typedef struct { pcg128_t state; } pcg_state_128;

typedef struct {
  pcg128_t state;
  pcg128_t inc;
} pcg_state_setseq_128;

#define PCG_DEFAULT_MULTIPLIER_HIGH 2549297995355413924ULL
#define PCG_DEFAULT_MULTIPLIER_LOW 4865540595714422341ULL

#define PCG_DEFAULT_MULTIPLIER_128                                             \
  PCG_128BIT_CONSTANT(PCG_DEFAULT_MULTIPLIER_HIGH, PCG_DEFAULT_MULTIPLIER_LOW)
#define PCG_DEFAULT_INCREMENT_128                                              \
  PCG_128BIT_CONSTANT(6364136223846793005ULL, 1442695040888963407ULL)
#define PCG_STATE_SETSEQ_128_INITIALIZER                                       \
  {                                                                            \
    PCG_128BIT_CONSTANT(0x979c9a98d8462005ULL, 0x7d3e9cb6cfe0549bULL)          \
    , PCG_128BIT_CONSTANT(0x0000000000000001ULL, 0xda3e39cb94b95bdbULL)        \
  }

static inline uint64_t pcg_rotr_64(uint64_t value, unsigned int rot) {
#ifdef _WIN32
  return _rotr64(value, rot);
#else
  return (value >> rot) | (value << ((-rot) & 63));
#endif
}

#ifdef PCG_EMULATED_128BIT_MATH

static inline pcg128_t pcg128_add(pcg128_t a, pcg128_t b) {
  pcg128_t result;

  result.low = a.low + b.low;
  result.high = a.high + b.high + (result.low < b.low);
  return result;
}

static inline void _pcg_mult64(uint64_t x, uint64_t y, uint64_t *z1,
                               uint64_t *z0) {

#if defined _WIN32 && _MSC_VER >= 1900 && _M_AMD64
  z0[0] = _umul128(x, y, z1);
#else
  uint64_t x0, x1, y0, y1;
  uint64_t w0, w1, w2, t;
  /* Lower 64 bits are straightforward clock-arithmetic. */
  *z0 = x * y;

  x0 = x & 0xFFFFFFFFULL;
  x1 = x >> 32;
  y0 = y & 0xFFFFFFFFULL;
  y1 = y >> 32;
  w0 = x0 * y0;
  t = x1 * y0 + (w0 >> 32);
  w1 = t & 0xFFFFFFFFULL;
  w2 = t >> 32;
  w1 += x0 * y1;
  *z1 = x1 * y1 + w2 + (w1 >> 32);
#endif
}

static inline pcg128_t pcg128_mult(pcg128_t a, pcg128_t b) {
  uint64_t h1;
  pcg128_t result;

  h1 = a.high * b.low + a.low * b.high;
  _pcg_mult64(a.low, b.low, &(result.high), &(result.low));
  result.high += h1;
  return result;
}

static inline void pcg_setseq_128_step_r(pcg_state_setseq_128 *rng) {
  rng->state = pcg128_add(pcg128_mult(rng->state, PCG_DEFAULT_MULTIPLIER_128),
                           rng->inc);
}

static inline uint64_t pcg_output_xsl_rr_128_64(pcg128_t state) {
  return pcg_rotr_64(state.high ^ state.low, state.high >> 58u);
}

static inline void pcg_setseq_128_srandom_r(pcg_state_setseq_128 *rng,
                                            pcg128_t initstate,
                                            pcg128_t initseq) {
  rng->state = PCG_128BIT_CONSTANT(0ULL, 0ULL);
  rng->inc.high = initseq.high << 1u;
  rng->inc.high |= initseq.low >> 63u;
  rng->inc.low = (initseq.low << 1u) | 1u;
  pcg_setseq_128_step_r(rng);
  rng->state = pcg128_add(rng->state, initstate);
  pcg_setseq_128_step_r(rng);
}

static inline uint64_t
pcg_setseq_128_xsl_rr_64_random_r(pcg_state_setseq_128 *rng) {
#if defined _WIN32 && _MSC_VER >= 1900 && _M_AMD64
  uint64_t h1;
  pcg128_t product;

  /* Manually inline the multiplication and addition using intrinsics */
  h1 = rng->state.high * PCG_DEFAULT_MULTIPLIER_LOW +
       rng->state.low * PCG_DEFAULT_MULTIPLIER_HIGH;
  product.low =
      _umul128(rng->state.low, PCG_DEFAULT_MULTIPLIER_LOW, &(product.high));
  product.high += h1;
  _addcarry_u64(_addcarry_u64(0, product.low, rng->inc.low, &(rng->state.low)),
                product.high, rng->inc.high, &(rng->state.high));
  return _rotr64(rng->state.high ^ rng->state.low, rng->state.high >> 58u);
#else
  pcg_setseq_128_step_r(rng);
  return pcg_output_xsl_rr_128_64(rng->state);
#endif
}

#else /* PCG_EMULATED_128BIT_MATH */

static inline void pcg_setseq_128_step_r(pcg_state_setseq_128 *rng) {
  rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_128 + rng->inc;
}

static inline uint64_t pcg_output_xsl_rr_128_64(pcg128_t state) {
  return pcg_rotr_64(((uint64_t)(state >> 64u)) ^ (uint64_t)state,
                     state >> 122u);
}

static inline uint64_t
pcg_setseq_128_xsl_rr_64_random_r(pcg_state_setseq_128* rng)
{
    pcg_setseq_128_step_r(rng);
    return pcg_output_xsl_rr_128_64(rng->state);
}

static inline void pcg_setseq_128_srandom_r(pcg_state_setseq_128 *rng,
                                            pcg128_t initstate,
                                            pcg128_t initseq) {
  rng->state = 0U;
  rng->inc = (initseq << 1u) | 1u;
  pcg_setseq_128_step_r(rng);
  rng->state += initstate;
  pcg_setseq_128_step_r(rng);
}

#endif /* PCG_EMULATED_128BIT_MATH */

static inline uint64_t
pcg_setseq_128_xsl_rr_64_boundedrand_r(pcg_state_setseq_128 *rng,
                                       uint64_t bound) {
  uint64_t threshold = -bound % bound;
  for (;;) {
    uint64_t r = pcg_setseq_128_xsl_rr_64_random_r(rng);
    if (r >= threshold)
      return r % bound;
  }
}

extern pcg128_t pcg_advance_lcg_128(pcg128_t state, pcg128_t delta,
                                    pcg128_t cur_mult, pcg128_t cur_plus);

static inline void pcg_setseq_128_advance_r(pcg_state_setseq_128 *rng,
                                            pcg128_t delta) {
  rng->state = pcg_advance_lcg_128(rng->state, delta,
                                   PCG_DEFAULT_MULTIPLIER_128, rng->inc);
}

typedef pcg_state_setseq_128 pcg64_random_t;
#define pcg64_random_r pcg_setseq_128_xsl_rr_64_random_r
#define pcg64_boundedrand_r pcg_setseq_128_xsl_rr_64_boundedrand_r
#define pcg64_srandom_r pcg_setseq_128_srandom_r
#define pcg64_advance_r pcg_setseq_128_advance_r
#define PCG64_INITIALIZER PCG_STATE_SETSEQ_128_INITIALIZER

#if __cplusplus
}
#endif

typedef struct s_pcg64_state {
  pcg64_random_t *pcg_state;
  int has_uint32;
  uint32_t uinteger;
} pcg64_state;

static inline uint64_t pcg64_next64(pcg64_state *state) {
  return pcg64_random_r(state->pcg_state);
}

static inline uint32_t pcg64_next32(pcg64_state *state) {
  uint64_t next;
  if (state->has_uint32) {
    state->has_uint32 = 0;
    return state->uinteger;
  }
  next = pcg64_random_r(state->pcg_state);
  state->has_uint32 = 1;
  state->uinteger = (uint32_t)(next >> 32);
  return (uint32_t)(next & 0xffffffff);
}

void pcg64_advance(pcg64_state *state, uint64_t *step);

void pcg64_set_seed(pcg64_state *state, uint64_t *seed, uint64_t *inc);

void pcg64_get_state(pcg64_state *state, uint64_t *state_arr, int *has_uint32,
                     uint32_t *uinteger);

void pcg64_set_state(pcg64_state *state, uint64_t *state_arr, int has_uint32,
                     uint32_t uinteger);

#endif /* PCG64_H_INCLUDED */

Read our documentation on viewing source code .

Navigation	Overlay
`t` Navigate files	`h` Toggle hits
`y` Change url to tip of branch	`m` Toggle misses
`b / v` Jump to prev/next hit line	`p` Toggle partial
`z / x` Jump to prev/next missed or partial line	`1..9` Toggle flags
`shift + o` Open current page in GitHub	`a` Toggle all on
`/ or ?` Show keyboard shortcuts dialog	`c` Toggle context lines or commits

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