diff options
| -rw-r--r-- | benchmarks/crand_benchmark2.cpp | 25 | ||||
| -rw-r--r-- | docs/clist_api.md | 1 | ||||
| -rw-r--r-- | docs/crand_api.md | 6 | ||||
| -rw-r--r-- | stc/crand.h | 17 |
4 files changed, 41 insertions, 8 deletions
diff --git a/benchmarks/crand_benchmark2.cpp b/benchmarks/crand_benchmark2.cpp index 966a3675..2692e760 100644 --- a/benchmarks/crand_benchmark2.cpp +++ b/benchmarks/crand_benchmark2.cpp @@ -4,6 +4,28 @@ #include "stc/crand.h"
#include "others/pcg_random.hpp"
+static struct stc32_state { stc64_t rng; uint64_t spare; unsigned n; } stc32_global =
+ {{0x7a5fed, 0x8e3f52, 0x9bc713, 0x6a09e667a7541669}, 0, 0};
+
+STC_INLINE void stc32_srandom(uint64_t seed) { stc32_global.rng = stc64_init(seed); }
+STC_INLINE uint32_t stc32_random(void) {
+ return (uint32_t) (++stc32_global.n & 1 ? (stc32_global.spare = stc64_rand(&stc32_global.rng))
+ : (stc32_global.spare >> 32));
+}
+
+static unsigned long myrand_next = 1;
+
+/* RAND_MAX assumed to be 32767 */
+int myrand(void) {
+ myrand_next = myrand_next * 214013 + 2531011;
+ return (myrand_next >> 16) & 0x7fff;
+}
+
+void mysrand(unsigned seed) {
+ myrand_next = seed;
+}
+
+
enum {N = 1000000000};
void test1(void)
@@ -88,7 +110,8 @@ void test3(void) before = clock();
sum = 0;
c_forrange (N) {
- sum += stc64_rand(&rng);
+ //sum += stc64_rand(&rng);
+ sum += rand();
}
diff = clock() - before;
printf("stc64_random:\t\t%.02f, %zu sz:%zu\n", (float) diff / CLOCKS_PER_SEC, sum, sizeof rng);
diff --git a/docs/clist_api.md b/docs/clist_api.md index deaec104..b056629d 100644 --- a/docs/clist_api.md +++ b/docs/clist_api.md @@ -95,6 +95,7 @@ clist_X_iter_t clist_X_begin(const clist_X* self); clist_X_iter_t clist_X_end(const clist_X* self); void clist_X_next(clist_X_iter_t* it); clist_X_value_t* clist_X_itval(clist_X_iter_t it); +clist_X_iter_t clist_X_fwd(clist_X_iter it, size_t n); clist_X_value_t clist_X_value_clone(clist_X_value_t val); ``` diff --git a/docs/crand_api.md b/docs/crand_api.md index 46d6ef1c..995c5fa0 100644 --- a/docs/crand_api.md +++ b/docs/crand_api.md @@ -40,14 +40,20 @@ All cstr definitions and prototypes may be included in your C source file by inc ## Methods ```c + void stc64_srandom(uint64_t seed); + uint64_t stc64_random(void); + 1) stc64_t stc64_init(uint64_t seed); 2) stc64_t stc64_with_seq(uint64_t seed, uint64_t seq); + 3) uint64_t stc64_rand(stc64_t* rng); 4) double stc64_randf(stc64_t* rng); + 5) stc64_uniform_t stc64_uniform_init(int64_t low, int64_t high); 6) int64_t stc64_uniform(stc64_t* rng, stc64_uniform_t* dist); 7) stc64_uniformf_t stc64_uniformf_init(double low, double high); 8) double stc64_uniformf(stc64_t* rng, stc64_uniformf_t* dist); + 9) stc64_normalf_t stc64_normalf_init(double mean, double stddev); 10) double stc64_normalf(stc64_t* rng, stc64_normalf_t* dist); ``` diff --git a/stc/crand.h b/stc/crand.h index 3d010eb2..929634ee 100644 --- a/stc/crand.h +++ b/stc/crand.h @@ -46,7 +46,7 @@ int main() { typedef struct {uint64_t state[4];} stc64_t;
typedef struct {int64_t lower; uint64_t range, threshold;} stc64_uniform_t;
typedef struct {double lower, range;} stc64_uniformf_t;
-typedef struct {double mean, stddev, next; bool has_next;} stc64_normalf_t;
+typedef struct {double mean, stddev, next; unsigned has_next;} stc64_normalf_t;
/* Stc64: random number generator, range [0, 2^64). PRNG copyright Tyge Løvset, NORCE Research, 2020 */
@@ -61,6 +61,11 @@ STC_INLINE uint64_t stc64_rand(stc64_t* rng) { return result;
}
+/* Global random() */
+static stc64_t stc64_global = {{0x26aa069ea2fb1a4d, 0x70c72c95cd592d04, 0x504f333d3aa0b359, 0x6a09e667a754166b}};
+STC_INLINE void stc64_srandom(uint64_t seed) { stc64_global = stc64_init(seed); }
+STC_INLINE uint64_t stc64_random(void) { return stc64_rand(&stc64_global); }
+
/* Float64 random number in range [low, high). */
STC_INLINE double stc64_randf(stc64_t* rng) {
union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (stc64_rand(rng) >> 12)};
@@ -100,7 +105,7 @@ STC_INLINE int64_t stc64_uniform(stc64_t* rng, stc64_uniform_t* d) { /* Normal distributed RNG, Float64. */
STC_INLINE stc64_normalf_t stc64_normalf_init(double mean, double stddev) {
- stc64_normalf_t dist = {mean, stddev, 0.0, false}; return dist;
+ stc64_normalf_t dist = {mean, stddev, 0.0, 0}; return dist;
}
STC_API double stc64_normalf(stc64_t* rng, stc64_normalf_t* dist);
@@ -118,7 +123,7 @@ STC_API double stc64_normalf(stc64_t* rng, stc64_normalf_t* dist); */
STC_DEF stc64_t stc64_init(uint64_t seed) {
- return stc64_with_seq(seed, 0x3504f333d3aa0b34);
+ return stc64_with_seq(seed, seed + 0x3504f333d3aa0b34);
}
STC_DEF stc64_t stc64_with_seq(uint64_t seed, uint64_t seq) {
stc64_t rng = {{seed, seed, seed, (seq << 1u) | 1u}};
@@ -136,17 +141,15 @@ STC_DEF stc64_uniform_t stc64_uniform_init(int64_t low, int64_t high) { /* Marsaglia polar method for gaussian/normal distribution. */
STC_DEF double stc64_normalf(stc64_t* rng, stc64_normalf_t* dist) {
double u1, u2, s, m;
- if (dist->has_next) {
- dist->has_next = false;
+ if (dist->has_next++ & 1)
return dist->next * dist->stddev + dist->mean;
- }
do {
u1 = 2.0 * stc64_randf(rng) - 1.0;
u2 = 2.0 * stc64_randf(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, dist->has_next = true;
+ dist->next = u2 * m;
return (u1 * m) * dist->stddev + dist->mean;
}
|