-
Notifications
You must be signed in to change notification settings - Fork 71
/
crand.h
161 lines (140 loc) · 5.57 KB
/
crand.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
/* MIT License
*
* Copyright (c) 2023 Tyge Løvset
*
* 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.
*/
#include "priv/linkage.h"
#ifndef CRAND_H_INCLUDED
#define CRAND_H_INCLUDED
#include "ccommon.h"
/*
// crand: Pseudo-random number generator
#include "stc/crand.h"
int main(void) {
uint64_t seed = 123456789;
crand_t rng = crand_init(seed);
crand_uniform_t dist1 = crand_uniform_init(1, 6);
crand_normal_t dist3 = crand_normal_init(1.0, 10.0);
uint64_t i = crand_u64(&rng);
int64_t iu = crand_uniform(&rng, &dist1);
double xn = crand_normal(&rng, &dist3);
}
*/
#include <string.h>
#include <math.h>
typedef struct crand { uint64_t state[5]; } crand_t;
typedef struct crand_uniform { int64_t lower; uint64_t range, threshold; } crand_uniform_t;
typedef struct crand_normal { double mean, stddev, next; int has_next; } crand_normal_t;
/* PRNG crand_t.
* Very fast PRNG suited for parallel usage with Weyl-sequence parameter.
* 320-bit state, 256 bit is mutable.
* Noticable faster than xoshiro and pcg, slighly slower than wyrand64 and
* Romu, but these have restricted capacity for larger parallel jobs or unknown minimum periods.
* crand_t supports 2^63 unique threads with a minimum 2^64 period lengths each.
* Passes all statistical tests, e.g PractRand and correlation tests, i.e. interleaved
* streams with one-bit diff state. Even the 16-bit version (LR=6, RS=5, LS=3) passes
* PractRand to multiple TB input.
*/
/* Global crand_t PRNGs */
STC_API void csrand(uint64_t seed);
STC_API uint64_t crand(void);
STC_API double crandf(void);
/* Init crand_t prng with a seed */
STC_API crand_t crand_init(uint64_t seed);
/* Unbiased bounded uniform distribution. range [low, high] */
STC_API crand_uniform_t crand_uniform_init(int64_t low, int64_t high);
STC_API int64_t crand_uniform(crand_t* rng, crand_uniform_t* dist);
/* Normal/gaussian distribution. */
STC_INLINE crand_normal_t crand_normal_init(double mean, double stddev)
{ crand_normal_t r = {mean, stddev, 0.0, 0}; return r; }
STC_API double crand_normal(crand_t* rng, crand_normal_t* dist);
/* Main crand_t prng */
STC_INLINE uint64_t crand_u64(crand_t* rng) {
uint64_t *s = rng->state;
const uint64_t result = (s[0] ^ (s[3] += s[4])) + s[1];
s[0] = s[1] ^ (s[1] >> 11);
s[1] = s[2] + (s[2] << 3);
s[2] = ((s[2] << 24) | (s[2] >> (64 - 24))) + result;
return result;
}
/* Float64 random number in range [0.0, 1.0). */
STC_INLINE double crand_f64(crand_t* rng) {
union {uint64_t i; double f;} u = {0x3FF0000000000000U | (crand_u64(rng) >> 12)};
return u.f - 1.0;
}
/* -------------------------- IMPLEMENTATION ------------------------- */
#if defined(i_implement) || defined(i_static)
/* Global random seed */
static crand_t crand_global = {{ // csrand(0)
0x9e3779b97f4a7c15, 0x6f68261b57e7a770,
0xe220a838bf5c9dde, 0x7c17d1800457b1ba, 0x1,
}};
STC_DEF void csrand(uint64_t seed)
{ crand_global = crand_init(seed); }
STC_DEF uint64_t crand(void)
{ return crand_u64(&crand_global); }
STC_DEF double crandf(void)
{ return crand_f64(&crand_global); }
STC_DEF crand_t crand_init(uint64_t seed) {
crand_t rng; uint64_t* s = rng.state;
s[0] = seed + 0x9e3779b97f4a7c15;
s[1] = (s[0] ^ (s[0] >> 30))*0xbf58476d1ce4e5b9;
s[2] = (s[1] ^ (s[1] >> 27))*0x94d049bb133111eb;
s[3] = s[0] ^ s[2] ^ (s[2] >> 31);
s[4] = (seed << 1) | 1;
return rng;
}
/* Init unbiased uniform uint RNG with bounds [low, high] */
STC_DEF crand_uniform_t crand_uniform_init(int64_t low, int64_t high) {
crand_uniform_t dist = {low, (uint64_t) (high - low + 1)};
dist.threshold = (uint64_t)(0 - dist.range) % dist.range;
return dist;
}
/* Int64 uniform distributed RNG, range [low, high]. */
STC_DEF int64_t crand_uniform(crand_t* rng, crand_uniform_t* d) {
uint64_t lo, hi;
#ifdef c_umul128
do { c_umul128(crand_u64(rng), d->range, &lo, &hi); } while (lo < d->threshold);
#else
do { lo = crand_u64(rng); hi = lo % d->range; } while (lo - hi > -d->range);
#endif
return d->lower + (int64_t)hi;
}
/* Normal distribution PRNG. Marsaglia polar method */
STC_DEF double crand_normal(crand_t* rng, crand_normal_t* dist) {
double u1, u2, s, m;
if (dist->has_next++ & 1)
return dist->next*dist->stddev + dist->mean;
do {
u1 = 2.0 * crand_f64(rng) - 1.0;
u2 = 2.0 * crand_f64(rng) - 1.0;
s = u1*u1 + u2*u2;
} while (s >= 1.0 || s == 0.0);
m = sqrt(-2.0 * log(s) / s);
dist->next = u2*m;
return (u1*m)*dist->stddev + dist->mean;
}
#endif
#endif
#undef i_opt
#undef i_static
#undef i_header
#undef i_implement
#undef i_import