P99
p99_rand.h
Go to the documentation of this file.
1 /* This may look like nonsense, but it really is -*- mode: C; coding: utf-8 -*- */
2 /* */
3 /* Except for parts copied from previous work and as explicitly stated below, */
4 /* the author and copyright holder for this work is */
5 /* (C) copyright 2012-2014 Jens Gustedt, INRIA, France */
6 /* */
7 /* This file is free software; it is part of the P99 project. */
8 /* */
9 /* Licensed under the Apache License, Version 2.0 (the "License"); */
10 /* you may not use this file except in compliance with the License. */
11 /* You may obtain a copy of the License at */
12 /* */
13 /* http://www.apache.org/licenses/LICENSE-2.0 */
14 /* */
15 /* Unless required by applicable law or agreed to in writing, software */
16 /* distributed under the License is distributed on an "AS IS" BASIS, */
17 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
18 /* See the License for the specific language governing permissions and */
19 /* limitations under the License. */
20 /* */
21 #ifndef P99_RAND_H
22 #define P99_RAND_H
23 #include "p99_threads.h"
24 #include "p99_new.h"
25 #include "p99_clib.h"
26 
38 /* From: George Marsaglia (geo@stat.fsu.edu) */
39 /* Subject: Re: good C random number generator */
40 /* Newsgroups: comp.lang.c */
41 /* Date: 2003-05-13 08:55:05 PST */
42 
43 /* Here is an example with k=5, period about 2^160, one of the fastest
44  long period RNGs, returns more than 120 million random 32-bit
45  integers/second (1.8MHz CPU), seems to pass all tests: */
46 
47 P99_CONSTANT(int, p00_seed160_len, 6);
48 typedef uint32_t p00_seed160[p00_seed160_len];
49 
50 p99_inline
51 uint32_t p00_xorshift(p00_seed160 * p00_s) {
52  /* Use the 6th element as the running index, and compute the
53  necessary derived indices mod 5 */
54  ++(*p00_s)[5];
55  register uint32_t const p00_0 = ((*p00_s)[5]) % 5;
56  register uint32_t const p00_4 = p00_0 == 0 ? 4 : p00_0-1;
57  register uint32_t const p00_2 = (p00_0 == 4) ? 1 : ((p00_0 == 3) ? 0 : p00_0+2);
58 
59  /* Load the part of the state variable that we need. */
60  register uint32_t const p00_x = (*p00_s)[p00_0];
61  register uint32_t const p00_y = (*p00_s)[p00_2];
62  register uint32_t const p00_v = (*p00_s)[p00_4];
63 
64  register uint32_t t = (p00_x^(p00_x>>7));
65  t = (p00_v^(p00_v<<6))^(t^(t<<13));
66 
67  /* Change just one place in the state variables. */
68  (*p00_s)[p00_4] = t;
69 
70  /* Multiply with an odd number */
71  return (2*p00_y + 1) * t;
72 }
73 
117 typedef p00_seed160 p99_seed[2];
118 
119 struct p00_rand160 {
120  p99_once_flag p00_flag;
121  p99_seed p00_seed;
122 };
123 
124 P99_DECLARE_THREAD_LOCAL(struct p00_rand160, p00_seed);
125 
126 #define P00_BIGPRIME \
127 UINT64_C(10007814641597694113), \
128 UINT64_C(10015183610531627897), \
129 UINT64_C(10089390291074425231), \
130 UINT64_C(10117275823191396191), \
131 UINT64_C(10215588060907623179), \
132 UINT64_C(10389441874873414061), \
133 UINT64_C(10602620510410479149), \
134 UINT64_C(10690793455755991027), \
135 UINT64_C(10793439684376201283), \
136 UINT64_C(10969167379420052431), \
137 UINT64_C(11050883422537956197), \
138 UINT64_C(11110712460003231287), \
139 UINT64_C(11342254921836171103), \
140 UINT64_C(11388752923558666351), \
141 UINT64_C(11682943258734137249), \
142 UINT64_C(11697622259988783581), \
143 UINT64_C(11731205802697935733), \
144 UINT64_C(11799524298339249581), \
145 UINT64_C(11862963432241722239), \
146 UINT64_C(11903172551239393097), \
147 UINT64_C(11937270107112816793), \
148 UINT64_C(11962673972948293321), \
149 UINT64_C(12001454206321043837), \
150 UINT64_C(12066567988703197129), \
151 UINT64_C(12090279337620046961), \
152 UINT64_C(12320220036580238077), \
153 UINT64_C(12409336420886496139), \
154 UINT64_C(12441355201341188273), \
155 UINT64_C(12444926949253327381), \
156 UINT64_C(12464345941324672183), \
157 UINT64_C(12469071720100075169), \
158 UINT64_C(12541090550165194183), \
159 UINT64_C(12558029109792848393), \
160 UINT64_C(12629045726855953273), \
161 UINT64_C(12845265184712635099), \
162 UINT64_C(12904165305535099733), \
163 UINT64_C(12949777725323667257), \
164 UINT64_C(13014679980528315931), \
165 UINT64_C(13027415158206533719), \
166 UINT64_C(13073256976558063279), \
167 UINT64_C(13153853421202095331), \
168 UINT64_C(13381351390611252151), \
169 UINT64_C(13825543164380321747), \
170 UINT64_C(13917403769616578911), \
171 UINT64_C(13929904570285153753), \
172 UINT64_C(14140216129016347127), \
173 UINT64_C(14184039880925065139), \
174 UINT64_C(14271570646513640671), \
175 UINT64_C(14446852824178510313), \
176 UINT64_C(14666997777683867003), \
177 UINT64_C(14692182605128454039), \
178 UINT64_C(14771388299183335229), \
179 UINT64_C(15059759696658839237), \
180 UINT64_C(15091316017205268437), \
181 UINT64_C(15180626055737690959), \
182 UINT64_C(15339437418814258573), \
183 UINT64_C(15601183704045111169), \
184 UINT64_C(15983170142123181797), \
185 UINT64_C(16545722488476290101), \
186 UINT64_C(16580005714519352107), \
187 UINT64_C(16589193612296178103), \
188 UINT64_C(16620171437104037921), \
189 UINT64_C(16713706803710883721), \
190 UINT64_C(16826711990149094791), \
191 UINT64_C(16886322755473635461), \
192 UINT64_C(16950140700365130619), \
193 UINT64_C(16952043704640877837), \
194 UINT64_C(17020386452131177189), \
195 UINT64_C(17020824347466355583), \
196 UINT64_C(17044414685455096133), \
197 UINT64_C(17089443963053420461), \
198 UINT64_C(17165479835541081871), \
199 UINT64_C(17185820514261086599), \
200 UINT64_C(17191892745505804273), \
201 UINT64_C(17613198595545138731), \
202 UINT64_C(17653178755564367203), \
203 UINT64_C(17653752776501147281), \
204 UINT64_C(17728444844615762171), \
205 UINT64_C(17744620187156425403), \
206 UINT64_C(17816722786105806973), \
207 UINT64_C(17831662963314755641), \
208 UINT64_C(17886004545482299117), \
209 UINT64_C(17920151933265509833), \
210 UINT64_C(17969552600607433963), \
211 UINT64_C(18066918898331824901), \
212 UINT64_C(18094173144238831753), \
213 UINT64_C(18134240252569387847), \
214 UINT64_C(18219113917191524677), \
215 UINT64_C(18298168206731166317), \
216 UINT64_C(18319234190200763803)
217 
218 P99_CONST_FUNCTION
219 p99_inline
220 uint32_t p00_bitpack(void const* p00_p) {
221  uintptr_t p00_u = (uintptr_t)p00_p;
222 #if UINTPTR_MAX == UINT32_MAX
223  return p00_u;
224 #else
225  uint32_t p00_r = 0;
226  do {
227  p00_r ^= p00_u;
228  p00_u >>= (sizeof(uint32_t)*CHAR_BIT);
229  } while (p00_u);
230  return p00_r;
231 #endif
232 }
233 
234 
235 P99_WEAK(p00_rand_init)
236 void p00_rand_init(void* p00_p) {
237  p99_seed * p00_s = p00_p
238  ? p00_p
239  : &P99_THREAD_LOCAL(p00_seed).p00_seed;
240  /* local to this initialization call, stack address */
241  uint32_t p00_0 = p00_bitpack(&p00_s);
242  /* unique to this thread, heap address */
243  uint32_t p00_1 = p00_bitpack(p00_p);
244  struct timespec p00_ts;
245  timespec_get(&p00_ts, TIME_UTC);
246  uint32_t p00_2 = p00_ts.tv_sec;
247  uint32_t p00_3 = p00_ts.tv_nsec;
248 #ifdef TIME_MONOTONIC
249  struct timespec p00_tm;
250  timespec_get(&p00_tm, TIME_MONOTONIC);
251  uint32_t p00_4 = p00_tm.tv_nsec;
252 #else
253  uint32_t p00_4 = p00_bitpack(__func__);
254 #endif
255  /* index, unique to this thread */
256  uint32_t p00_ind = p00_0 ^ p00_1 ^ p00_2 ^ p00_3 ^ p00_4;
257  p00_seed160 p00_st = {
258  p00_0,
259  p00_1,
260  p00_2,
261  p00_3,
262  p00_4,
263  p00_ind,
264  };
265  /* mix things up a bit */
266  for (unsigned p00_i = 0; p00_i < 32; ++p00_i) p00_xorshift(&p00_st);
267  /* Now create two different state vectors that are de-phased. We
268  can be in any state of the xorshift generator, so there are
269  2^160 different initializations. */
270  for (unsigned p00_j = 0; p00_j < 2; ++p00_j) {
271  for (unsigned p00_i = 0; p00_i < p00_seed160_len; ++p00_i)
272  (*p00_s)[p00_j][p00_i] = p00_xorshift(&p00_st);
273  }
274 }
275 
308 p99_inline
309 p99_seed * p99_seed_get(void) {
310  struct p00_rand160 * p00_loc = &P99_THREAD_LOCAL(p00_seed);
311  p99_call_once(&p00_loc->p00_flag, p00_rand_init, &p00_loc->p00_seed);
312  return &p00_loc->p00_seed;
313 }
314 
315 P99_WEAK(p00_bigprime)
316 uint64_t const p00_bigprime[] = { P00_BIGPRIME };
317 
318 P99_CONSTANT(int, p00_bigprime_len, P99_ALEN(p00_bigprime));
319 
338 P99_DEFARG_DOCU(p99_rand)
339 p99_inline
340 uint64_t p99_rand(register p99_seed * p00_seed) {
341  uint32_t p00_0 = p00_xorshift(&(*p00_seed)[0]);
342  uint64_t p00_1 = p00_xorshift(&(*p00_seed)[1]);
343  uint64_t p00_0r = p00_0 % p00_bigprime_len;
344  uint64_t p00_0d = p00_0 / p00_bigprime_len;
345  /* Use part of the bits to choose a big number */
346  p00_1 *= p00_bigprime[p00_0r];
347  /* Use the rest of the bits to add more randomness */
348  return (p00_0d ^ p00_1);
349 }
350 
351 #ifndef DOXYGEN
352 #define p99_rand(...) P99_CALL_DEFARG(p99_rand, 1, __VA_ARGS__)
353 #define p99_rand_defarg_0() (p99_seed_get())
354 #endif
355 
370 P99_DEFARG_DOCU(p99_drand)
371 p99_inline
372 double p99_drand(register p99_seed * p00_seed) {
373  enum {
374  p00_s = (DBL_MANT_DIG > 64 ? 64 : DBL_MANT_DIG),
375  p00_0 = 64 - p00_s,
376  p00_m = (1 << p00_0),
377  };
378  double const p00_imax = 1.0 / (UINT64_C(1) << p00_s);
379  uint64_t p00_r = p99_rand(p00_seed);
380  double p00_1 = ((p00_r >> p00_0)^(p00_r & p00_m)) * p00_imax;
381  return p00_1;
382 }
383 
384 #ifndef DOXYGEN
385 #define p99_drand(...) P99_CALL_DEFARG(p99_drand, 1, __VA_ARGS__)
386 #define p99_drand_defarg_0() (p99_seed_get())
387 #endif
388 
393 #endif
TIME_UTC
#define TIME_UTC
expands to an integer constant greater than 0 that designates the UTC time base since an implementati...
Definition: p99_clib.h:197
p99_new.h
Macros for initialization and allocation.
p99_once_flag
complete object type that holds a flag for use by p99_call_once
Definition: p99_threads.h:62
P99_THREAD_LOCAL
#define P99_THREAD_LOCAL
P99_CONSTANT
#define P99_CONSTANT(T, NAME, INIT)
define a compile time constant NAME of type T with value INIT
Definition: p99_enum.h:258
TIME_MONOTONIC
#define TIME_MONOTONIC
expands to an integer constant greater than 0 that designates a real time clock who's base is usually...
Definition: p99_clib.h:215
P99_ALEN
#define P99_ALEN(ARR, N)
Produce the length of the argument array ARR in terms of number of elements.
Definition: p99_for.h:678
p99_inline
#define p99_inline
Try to force a function always to be inlined.
Definition: p99_compiler.h:496
P99_WEAK
#define P99_WEAK(...)
Declare a symbol to be weak such that it can be provided several times without error.
Definition: p99_compiler.h:561
P99_DEFARG_DOCU
#define P99_DEFARG_DOCU(NAME)
Provide a documentation section to a function defined with P99_CALL_DEFARG.
Definition: p99_defarg.h:318
timespec_get
int timespec_get(struct timespec *p00_ts, int p00_base)
The timespec_get function sets the interval pointed to by p00_ts to hold the current calendar time ba...
Definition: p99_clib.h:260
P99_DECLARE_THREAD_LOCAL
#define P99_DECLARE_THREAD_LOCAL
P99_CONST_FUNCTION
#define P99_CONST_FUNCTION
On architectures that support this, assert that a function is "const", i.e only depends on parameters...
Definition: p99_compiler.h:622
p99_seed
p00_seed160 p99_seed[2]
The internal state that the PRG for ::p99_rand and ::p99_drand uses.
Definition: p99_rand.h:117
p99_threads.h
p99_clib.h
p99_call_once
#define p99_call_once(FLAG, FUNC, ARG)
Call a function FUNC exactly once, optionally providing it with argument ARG.
Definition: p99_threads.h:202
Generated by   doxygen 1.8.17