diff options
| -rw-r--r-- | benchmarks/cmap_benchmark.cpp (renamed from examples/benchmark.cpp) | 0 | ||||
| -rw-r--r-- | benchmarks/cpque_benchmark.cpp (renamed from examples/heap.c) | 2 | ||||
| -rw-r--r-- | benchmarks/crand_benchmark.cpp | 155 | ||||
| -rw-r--r-- | benchmarks/crand_benchmark2.cpp | 75 | ||||
| -rw-r--r-- | examples/birthday.c | 46 | ||||
| -rw-r--r-- | examples/rngtest.c | 40 | ||||
| -rw-r--r-- | stc/cdeq.h | 2 | ||||
| -rw-r--r-- | stc/crand.h | 31 |
8 files changed, 264 insertions, 87 deletions
diff --git a/examples/benchmark.cpp b/benchmarks/cmap_benchmark.cpp index a159d24b..a159d24b 100644 --- a/examples/benchmark.cpp +++ b/benchmarks/cmap_benchmark.cpp diff --git a/examples/heap.c b/benchmarks/cpque_benchmark.cpp index c7292d2e..625bb056 100644 --- a/examples/heap.c +++ b/benchmarks/cpque_benchmark.cpp @@ -11,7 +11,7 @@ int main() {
uint32_t seed = time(NULL);
crand_t rng;
- int N = 3000000, M = 100;
+ int N = 10000000, M = 10;
cpque_f pq = cpque_f_init();
diff --git a/benchmarks/crand_benchmark.cpp b/benchmarks/crand_benchmark.cpp new file mode 100644 index 00000000..78099905 --- /dev/null +++ b/benchmarks/crand_benchmark.cpp @@ -0,0 +1,155 @@ +#include <cstdint>
+#include <iostream>
+#include <ctime>
+#include <stc/crand.h>
+
+static inline uint64_t rotl64(const uint64_t x, const int k)
+ { return (x << k) | (x >> (64 - k)); }
+
+static uint64_t splitmix64_x = 87213627321ull; /* The state can be seeded with any value. */
+
+uint64_t splitmix64(void) {
+ uint64_t z = (splitmix64_x += 0x9e3779b97f4a7c15);
+ z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+ z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+ return z ^ (z >> 31);
+}
+
+static void init_state(uint64_t *rng, uint64_t seed) {
+ splitmix64_x = seed;
+ for (int i=0; i<4; ++i) rng[i] = splitmix64();
+}
+
+/* jsf64 */
+
+static inline uint64_t jsf64(uint64_t *s) {
+ uint64_t e = s[0] - rotl64(s[1], 7);
+ s[0] = s[1] ^ rotl64(s[2], 13);
+ s[1] = s[2] + rotl64(s[3], 37);
+ s[2] = s[3] + e;
+ s[3] = e + s[0];
+ return s[3];
+}
+
+/* xoshiro256** */
+
+static inline uint64_t xoshiro256starstar(uint64_t* s) {
+ const uint64_t result = rotl64(s[1] * 5, 7) * 9;
+ const uint64_t t = s[1] << 17;
+ s[2] ^= s[0];
+ s[3] ^= s[1];
+ s[1] ^= s[2];
+ s[0] ^= s[3];
+ s[2] ^= t;
+ s[3] = rotl64(s[3], 45);
+ return result;
+}
+
+
+inline unsigned long long wyhash64(uint64_t* s) {
+ *s += 0x60bee2bee120fc15ull;
+ __uint128_t tmp = (__uint128_t) *s * 0xa3b195354a39b70dull;
+ unsigned long long m1 = (tmp >> 64) ^ tmp;
+ tmp = (__uint128_t)m1 * 0x1b03738712fad5c9ull;
+ return (tmp >> 64) ^ tmp;
+}
+
+inline unsigned long long lehmer64(uint64_t* s) {
+ *(__uint128_t *)s *= 0xda942042e4dd58b5ull;
+ return *(__uint128_t *)s >> 64;
+}
+
+using namespace std;
+
+int main(void)
+{
+ enum {N = 524288000};
+ uint64_t* recipient = new uint64_t[N];
+ static crand_t rng;
+ init_state(rng.state, 12345123);
+
+ clock_t beg, end;
+ for (size_t ti = 0; ti < 4; ti++) {
+ beg = clock();
+ for (size_t i = 0; i < N; i++)
+ recipient[i] = wyhash64(rng.state);
+ end = clock();
+ cerr << "ROUND " << ti+1 << endl
+ << "wyhash64:\t"
+ << (float(end - beg) / CLOCKS_PER_SEC)
+ << " s" << endl;
+ cout << "bogus:" << recipient[312] << endl;
+ beg = clock();
+ for (size_t i = 0; i < N; i++)
+ recipient[i] = crand_next(&rng);
+ end = clock();
+ cerr << "stc crand:\t"
+ << (float(end - beg) / CLOCKS_PER_SEC)
+ << " s" << endl;
+ cout << "bogus:" << recipient[312] << endl;
+
+ beg = clock();
+ for (size_t i = 0; i < N; i++)
+ recipient[i] = xoshiro256starstar(rng.state);
+ end = clock();
+ cerr << "xoshiro256**:\t"
+ << (float(end - beg) / CLOCKS_PER_SEC)
+ << " s" << endl;
+ cout << "bogus:" << recipient[312] << endl;
+
+ beg = clock();
+ for (size_t i = 0; i < N; i++)
+ recipient[i] = lehmer64(rng.state);
+ end = clock();
+ cerr << "lehmer64:\t"
+ << ((float) end - beg) / CLOCKS_PER_SEC
+ << " s" << endl;
+ cout << "bogus:" << recipient[312] << endl;
+
+
+
+ cout << endl
+ << "Next we do random number computations only, doing no work."
+ << endl;
+ uint64_t s = 0;
+ beg = clock();
+ for (size_t i = 0; i < N; i++)
+ s += wyhash64(rng.state);
+ end = clock();
+ cerr << "wyhash64:\t"
+ << ((float) end - beg) / CLOCKS_PER_SEC
+ << " s" << endl;
+ cout << "bogus:" << s << endl;
+
+ s = 0;
+ beg = clock();
+ for (size_t i = 0; i < N; i++)
+ s += crand_next(&rng);
+ end = clock();
+ cerr << "stc crand:\t"
+ << ((float) end - beg) / CLOCKS_PER_SEC
+ << " s" << endl;
+ cout << "bogus:" << s << endl;
+
+ s = 0;
+ beg = clock();
+ for (size_t i = 0; i < N; i++)
+ s += xoshiro256starstar(rng.state);
+ end = clock();
+ cerr << "xoshiro256**:\t"
+ << ((float) end - beg) / CLOCKS_PER_SEC
+ << " s" << endl;
+ cout << "bogus:" << s << endl;
+
+ beg = clock();
+ for (size_t i = 0; i < N; i++)
+ s += lehmer64(rng.state);
+ end = clock();
+ cerr << "lehmer64:\t"
+ << ((float) end - beg) / CLOCKS_PER_SEC
+ << " s" << endl;
+ cout << "bogus:" << s << endl << endl;
+ }
+ delete[] recipient;
+ return 0;
+}
diff --git a/benchmarks/crand_benchmark2.cpp b/benchmarks/crand_benchmark2.cpp new file mode 100644 index 00000000..3ae6f8ab --- /dev/null +++ b/benchmarks/crand_benchmark2.cpp @@ -0,0 +1,75 @@ +#include <stdio.h>
+#include <time.h>
+#include <stc/crand.h>
+#include <random>
+
+void test1(void)
+{
+ enum {N = 1000000000};
+ clock_t diff, before;
+ uint64_t sum;
+
+ std::random_device device;
+ std::mt19937 rng(device());
+ std::uniform_int_distribution<int> idist(1, 10);
+ std::uniform_real_distribution<double> fdist(1, 10);
+
+ before = clock();
+ sum = 0;
+ c_forrange (N) {
+ sum += rng();
+ }
+ diff = clock() - before;
+ printf("std::random:\t\t%.02f, %zu\n", (float) diff / CLOCKS_PER_SEC, sum);
+
+ before = clock();
+ sum = 0;
+ c_forrange (N) {
+ sum += idist(rng);
+ }
+ diff = clock() - before;
+ printf("std::uniform:\t\t%.02f, %zu\n\n", (float) diff / CLOCKS_PER_SEC, sum);
+
+ c_forrange (30) printf("%02zd ", idist(rng));
+ puts("");
+ c_forrange (8) printf("%f ", fdist(rng));
+ puts("\n");
+}
+
+void test2(void)
+{
+ enum {N = 1000000000};
+ clock_t diff, before;
+ uint64_t sum;
+
+ crand_t rng = crand_init(time(NULL));
+ crand_uniform_t idist = crand_uniform_init(1, 10);
+ crand_uniformf_t fdist = crand_uniformf_init(1, 10);
+
+ before = clock();
+ sum = 0;
+ c_forrange (N) {
+ sum += crand_next(&rng);
+ }
+ diff = clock() - before;
+ printf("crand_next:\t\t%.02f, %zu\n", (float) diff / CLOCKS_PER_SEC, sum);
+
+ before = clock();
+ sum = 0;
+ c_forrange (N) {
+ sum += crand_uniform(&rng, &idist);
+ }
+ diff = clock() - before;
+ printf("crand_uniform:\t\t%.02f, %zu\n\n", (float) diff / CLOCKS_PER_SEC, sum);
+
+ c_forrange (30) printf("%02zd ", crand_uniform(&rng, &idist));
+ puts("");
+ c_forrange (8) printf("%f ", crand_uniformf(&rng, &fdist));
+ puts("\n");
+}
+
+int main()
+{
+ test1();
+ test2();
+}
\ No newline at end of file diff --git a/examples/birthday.c b/examples/birthday.c index 7ed30bf1..a2856a3f 100644 --- a/examples/birthday.c +++ b/examples/birthday.c @@ -1,61 +1,57 @@ -
+#include <math.h>
#include <stdio.h>
#include <time.h>
#include <stc/crand.h>
#include <stc/cmap.h>
-#include <stc/cvec.h>
using_cmap(ic, uint64_t, uint8_t);
+static uint64_t seed = 12345;
-const static uint64_t seed = 1234;
-const static uint64_t N = 1ull << 27;
-const static uint64_t mask = (1ull << 52) - 1;
-
-void repeats(void)
+static void test_repeats(void)
{
+ enum {BITS = 46, BITS_TEST = BITS/2 + 2};
+ const static uint64_t N = 1ull << BITS_TEST;
+ const static uint64_t mask = (1ull << BITS) - 1;
+
+ printf("birthday paradox: value range: 2^%d, testing repeats of 2^%d values\n", BITS, BITS_TEST);
crand_t rng = crand_init(seed);
cmap_ic m = cmap_ic_init();
cmap_ic_reserve(&m, N);
- clock_t now = clock();
c_forrange (i, N) {
uint64_t k = crand_next(&rng) & mask;
int v = ++cmap_ic_emplace(&m, k, 0).first->second;
- if (v > 1) printf("%zu: %llx - %d\n", i, k, v);
+ if (v > 1) printf("repeated value %llx (%d) at 2^%d\n", k, v, (int) log2(i));
}
- float diff = (float) (clock() - now) / CLOCKS_PER_SEC;
- printf("%.02f", diff);
}
using_cmap(x, uint32_t, uint64_t);
-using_cvec(x, uint64_t);
-void distribution(void)
+void test_distribution(void)
{
- crand_t rng = crand_init(seed); // time(NULL), time(NULL));
- const size_t N = 1ull << 28, M = 1ull << 9; // 1ull << 10;
- cmap_x map = cmap_x_with_capacity(M);
- clock_t now = clock();
- crand_uniform_t dist = crand_uniform_init(0, M);
+ enum {BITS = 26};
+ printf("distribution test: 2^%d values\n", BITS);
+ crand_t rng = crand_init(seed);
+ const size_t N = 1ull << BITS ;
+ cmap_x map = cmap_x_init();
c_forrange (N) {
- ++cmap_x_emplace(&map, crand_uniform(&rng, &dist), 0).first->second;
+ uint64_t k = crand_next(&rng);
+ ++cmap_x_emplace(&map, k & 0xf, 0).first->second;
}
- float diff = (float) (clock() - now) / CLOCKS_PER_SEC;
uint64_t sum = 0;
c_foreach (i, cmap_x, map) sum += i.ref->second;
sum /= map.size;
c_foreach (i, cmap_x, map)
- printf("%u: %zu - %zu\n", i.ref->first, i.ref->second, sum);
-
- printf("%.02f\n", diff);
+ printf("%4u: %zu - %zu: %11.8f\n", i.ref->first, i.ref->second, sum, (1 - (double) i.ref->second / sum));
}
int main()
{
- repeats();
- //distribution();
+ seed = time(NULL);
+ test_distribution();
+ test_repeats();
}
\ No newline at end of file diff --git a/examples/rngtest.c b/examples/rngtest.c deleted file mode 100644 index f52099c7..00000000 --- a/examples/rngtest.c +++ /dev/null @@ -1,40 +0,0 @@ -#include <stdio.h>
-#include <time.h>
-#include <stc/crand.h>
-#ifdef __cplusplus
-#include <random>
-#endif
-
-
-#define NN 1000000000
-
-int main(void)
-{
- clock_t diff, before;
- uint64_t v;
-
- crand_t rng = crand_init(time(NULL));
- crand_uniform_t idist = crand_uniform_init(10, 20);
- crand_uniformf_t fdist = crand_uniformf_init(10, 20);
-
- before = clock();
- v = 0;
- c_forrange (NN) {
- v += crand_next(&rng);
- }
- diff = clock() - before;
- printf("stc64_rand: %.02f, %zu\n", (float) diff / CLOCKS_PER_SEC, v);
-
- before = clock();
- v = 0;
- c_forrange (NN) {
- v += crand_uniform(&rng, &idist);
- }
- diff = clock() - before;
- printf("stc64_uniform: %.02f, %zu\n\n", (float) diff / CLOCKS_PER_SEC, v);
-
- c_forrange (30) printf("%02zd ", crand_uniform(&rng, &idist));
- puts("");
- c_forrange (8) printf("%f ", crand_uniformf(&rng, &fdist));
- puts("");
-}
\ No newline at end of file @@ -242,7 +242,7 @@ rep[0] = len, rep[1] = cap; \
self->base = (cdeq_##X##_value_t *) (rep + 2); \
self->data = self->base + nfront; \
- return _deq_##X##_expand(self, n, at_front); \
+ return _cdeq_##X##_expand(self, n, at_front); \
} \
size_t unused = cap - (len + n); \
size_t pos = at_front ? c_maxf(unused*0.9, (float) unused - nback) + n \
diff --git a/stc/crand.h b/stc/crand.h index 3652851c..ad477ea4 100644 --- a/stc/crand.h +++ b/stc/crand.h @@ -52,7 +52,15 @@ typedef struct {double mean, stddev, next; bool has_next;} crand_normalf_t; /* int random number generator, range [0, 2^64). PRNG copyright Tyge Løvset, NORCE Research, 2020 */
STC_API crand_t crand_init(uint64_t seed);
STC_API crand_t crand_with_seq(uint64_t seed, uint64_t seq);
-STC_API uint64_t crand_next(crand_t* rng);
+
+STC_INLINE uint64_t crand_next(crand_t* rng) {
+ enum {LROT = 24, RSHIFT = 11, LSHIFT = 3};
+ uint64_t *s = rng->state;
+ const uint64_t b = s[1], result = s[0] ^ (s[2] += s[3]|1);
+ s[0] = (b + (b << LSHIFT)) ^ (b >> RSHIFT);
+ s[1] = ((b << LROT) | (b >> (64 - LROT))) + result;
+ return result;
+}
/* double random number in range [low, high). */
STC_INLINE double crand_nextf(crand_t* rng) {
@@ -62,7 +70,6 @@ STC_INLINE double crand_nextf(crand_t* rng) { /* integer uniform distributed RNG, 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);
/* double uniform distributed RNG, range [low, high). */
STC_INLINE crand_uniformf_t crand_uniformf_init(double low, double high) {
@@ -85,9 +92,9 @@ STC_INLINE double crand_uniformf(crand_t* rng, crand_uniformf_t* dist) { : [lhs] "0" (a), [rhs] "rm" (b))
#endif
-STC_INLINE int64_t crand_uniform_fast(crand_t* rng, crand_uniform_t* d) {
+STC_INLINE int64_t crand_uniform(crand_t* rng, crand_uniform_t* d) {
uint64_t lo, hi;
- cmul128(crand_next(rng), d->range, &lo, &hi);
+ do { cmul128(crand_next(rng), d->range, &lo, &hi); } while (lo < d->threshold);
return d->lower + hi;
}
@@ -118,14 +125,6 @@ STC_DEF crand_t crand_with_seq(uint64_t seed, uint64_t seq) { for (int i = 0; i < 8; ++i) crand_next(&rng);
return rng;
}
-STC_DEF uint64_t crand_next(crand_t* rng) {
- enum {LROT = 24, RSHIFT = 11, LSHIFT = 3};
- uint64_t *s = rng->state;
- const uint64_t b = s[1], result = s[0] ^ (s[2] += s[3]|1);
- s[0] = (b + (b << LSHIFT)) ^ (b >> RSHIFT);
- s[1] = ((b << LROT) | (b >> (64 - LROT))) + result;
- return result;
-}
/* Very fast unbiased uniform int RNG with bounds [low, high] */
STC_DEF crand_uniform_t crand_uniform_init(int64_t low, int64_t high) {
@@ -134,14 +133,6 @@ STC_DEF crand_uniform_t crand_uniform_init(int64_t low, int64_t high) { return dist;
}
-STC_DEF int64_t crand_uniform(crand_t* rng, crand_uniform_t* d) {
- uint64_t lo, hi;
- do {
- cmul128(crand_next(rng), d->range, &lo, &hi);
- } while (lo < d->threshold);
- return d->lower + hi;
-}
-
/* Marsaglia polar method for gaussian/normal distribution. */
STC_DEF double crand_normalf(crand_t* rng, crand_normalf_t* dist) {
double u1, u2, s, m;
|