From 4b14954ce30ea7c5e63a49c50294bcec5978bbcd Mon Sep 17 00:00:00 2001 From: Tyge Løvset Date: Tue, 28 Jul 2020 14:58:20 +0200 Subject: Updated random --- examples/benchmark.c | 6 +- examples/heap.c | 8 +- examples/list.c | 4 +- examples/priority.c | 4 +- examples/rngbirthday.c | 8 +- examples/rngtest.c | 12 +-- stc/coptget.h | 173 ---------------------------------------- stc/crandom.h | 208 ++++++++++++------------------------------------- 8 files changed, 70 insertions(+), 353 deletions(-) delete mode 100644 stc/coptget.h diff --git a/examples/benchmark.c b/examples/benchmark.c index c9a82d6d..1c69b8ab 100644 --- a/examples/benchmark.c +++ b/examples/benchmark.c @@ -21,9 +21,9 @@ KHASH_MAP_INIT_INT64(ii, uint64_t) size_t seed; static const float maxLoadFactor = 0.77f; -crandom64_t rng; -#define SEED(s) rng = crandom64_init(seed) -#define RAND(N) (crandom64(&rng) & ((1 << N) - 1)) +CRand64 rng; +#define SEED(s) rng = crand64_init(seed) +#define RAND(N) (crand64_gen(&rng) & ((1 << N) - 1)) #define CMAP_SETUP(tag, Key, Value) CMap_##tag map = cmap_init \ diff --git a/examples/heap.c b/examples/heap.c index e4e91544..f0d56a46 100644 --- a/examples/heap.c +++ b/examples/heap.c @@ -9,12 +9,12 @@ declare_CVec_priority_queue(f, >); int main() { uint32_t seed = time(NULL); - crandom32_t pcg = crandom32_init(seed); + CRand32 pcg = crand32_init(seed); int N = 30000000, M = 100; CVec_f vec = cvec_init; clock_t start = clock(); for (int i=0; ivalue); else break; diff --git a/examples/priority.c b/examples/priority.c index c1890f8c..baab2779 100644 --- a/examples/priority.c +++ b/examples/priority.c @@ -9,12 +9,12 @@ declare_CVec(i, uint32_t); declare_CVec_priority_queue(i, >); // min-heap (increasing values) int main() { - crandom32_t pcg = crandom32_init(time(NULL)); + CRand32 pcg = crand32_init(time(NULL)); CVec_i heap = cvec_init; // Push ten million random numbers to queue for (int i=0; i<10000000; ++i) - cvecpq_i_push(&heap, crandom32(&pcg)); + cvecpq_i_push(&heap, crand32_gen(&pcg)); // Extract the hundred smallest. for (int i=0; i<100; ++i) { diff --git a/examples/rngbirthday.c b/examples/rngbirthday.c index 1e6178aa..8129497b 100644 --- a/examples/rngbirthday.c +++ b/examples/rngbirthday.c @@ -15,12 +15,12 @@ const static uint64_t mask = (1ull << 52) - 1; void repeats(void) { - crandom64_t rng = crandom64_init(seed); + CRand64 rng = crand64_init(seed); CMap_ic m = cmap_init; cmap_ic_reserve(&m, N); clock_t now = clock(); for (size_t i = 0; i < N; ++i) { - uint64_t k = crandom64(&rng) & mask; + uint64_t k = crand64_gen(&rng) & mask; int v = ++cmap_ic_insert(&m, k, 0)->value; if (v > 1) printf("%zu: %x - %d\n", i, k, v); } @@ -34,12 +34,12 @@ declare_CVec(x, uint64_t); void distribution(void) { - crandom32_t rng = crandom32_init(seed); // time(NULL), time(NULL)); + CRand32 rng = crand32_init(seed); // time(NULL), time(NULL)); const size_t N = 1ull << 28, M = 1ull << 9; // 1ull << 10; CMap_x map = cmap_x_make(M); clock_t now = clock(); for (size_t i = 0; i < N; ++i) { - ++cmap_x_insert(&map, crandom32b(&rng, M), 0)->value; + ++cmap_x_insert(&map, crand32_genBounded(&rng, M), 0)->value; } float diff = (float) (clock() - now) / CLOCKS_PER_SEC; diff --git a/examples/rngtest.c b/examples/rngtest.c index af7c7dca..a15a8fcb 100644 --- a/examples/rngtest.c +++ b/examples/rngtest.c @@ -14,26 +14,26 @@ int main(void) uint64_t v; printf("start\n"); - crandom32_t pcg = crandom32_init(time(NULL)); + CRand32 pcg = crand32_init(time(NULL)); before = clock(); \ v = 0; for (size_t i=0; ifaulty output field, and has a more consistent API. - - copt_get() is similar to GNU's getopt_long(). Each call parses one option and - returns the option name. opt->arg points to the option argument if present. - The function returns -1 when all command-line arguments are parsed. In this case, - opt->ind is the index of the first non-option argument. -Example: - int main(int argc, char *argv[]) - { - copt_long_t longopts[] = { - {"foo", copt_no_argument, 'f'}, - {"bar", copt_required_argument, 'b'}, - {"opt", copt_optional_argument, 'o'}, - {NULL} - }; - const char* optstr = "xy:z::123"; - printf("program -x -y ARG -z [ARG] -1 -2 -3 --foo --bar ARG --opt [ARG] [ARGUMENTS]\n"); - int c; - copt_t opt = copt_init; - while ((c = copt_get(&opt, argc, argv, optstr, longopts)) != -1) { - switch (c) { - case '?': printf("error: unknown option: %s\n", opt.faulty); break; - case ':': printf("error: missing argument for %s\n", opt.faulty); break; - default: printf("option: %c [%s]\n", c, opt.arg ? opt.arg : ""); break; - } - } - printf("\nNon-option arguments:"); - for (int i = opt.ind; i < argc; ++i) - printf(" %s", argv[i]); - putchar('\n'); - return 0; - } -*/ - -#include -#include - -enum { - copt_no_argument = 0, - copt_required_argument = 1, - copt_optional_argument = 2 -}; -typedef struct { - int ind; /* equivalent to optind */ - int opt; /* equivalent to optopt */ - char *arg; /* equivalent to optarg */ - char *faulty; /* points to the faulty option */ - int longindex; /* idx of long option; or -1 if short */ - int _i, _pos, _nargs; - char _faulty[4]; -} copt_t; - -typedef struct { - char *name; - int has_arg; - int val; -} copt_long_t; - -static const copt_t copt_init = {1, 0, NULL, NULL, -1, 1, 0, 0, {'-', '?', '\0'}}; - -static void _copt_permute(char *argv[], int j, int n) { /* move argv[j] over n elements to the left */ - int k; - char *p = argv[j]; - for (k = 0; k < n; ++k) - argv[j - k] = argv[j - k - 1]; - argv[j - k] = p; -} - -/* @param opt output; must be initialized to copt_init on first call - * @return ASCII val for a short option; longopt.val for a long option; - * -1 if argv[] is fully processed; '?' for an unknown option or - * an ambiguous long option; ':' if an option argument is missing - */ -static int copt_get(copt_t *opt, int argc, char *argv[], - const char *shortopts, const copt_long_t *longopts) { - int optc = -1, i0, j, posixly_correct = (shortopts[0] == '+'); - if (!posixly_correct) { - while (opt->_i < argc && (argv[opt->_i][0] != '-' || argv[opt->_i][1] == '\0')) - ++opt->_i, ++opt->_nargs; - } - opt->arg = 0, opt->longindex = -1, i0 = opt->_i; - if (opt->_i >= argc || argv[opt->_i][0] != '-' || argv[opt->_i][1] == '\0') { - opt->ind = opt->_i - opt->_nargs; - return -1; - } - if (argv[opt->_i][0] == '-' && argv[opt->_i][1] == '-') { /* "--" or a long option */ - if (argv[opt->_i][2] == '\0') { /* a bare "--" */ - _copt_permute(argv, opt->_i, opt->_nargs); - ++opt->_i, opt->ind = opt->_i - opt->_nargs; - return -1; - } - opt->opt = 0, optc = '?', opt->_pos = -1; - if (longopts) { /* parse long options */ - int k, n_exact = 0, n_partial = 0; - const copt_long_t *o = 0, *o_exact = 0, *o_partial = 0; - for (j = 2; argv[opt->_i][j] != '\0' && argv[opt->_i][j] != '='; ++j) {} /* find the end of the option name */ - for (k = 0; longopts[k].name != 0; ++k) - if (strncmp(&argv[opt->_i][2], longopts[k].name, j - 2) == 0) { - if (longopts[k].name[j - 2] == 0) ++n_exact, o_exact = &longopts[k]; - else ++n_partial, o_partial = &longopts[k]; - } - opt->faulty = argv[opt->_i]; - if (n_exact > 1 || (n_exact == 0 && n_partial > 1)) return '?'; - o = n_exact == 1? o_exact : n_partial == 1? o_partial : 0; - if (o) { - opt->opt = optc = o->val, opt->longindex = o - longopts; - if (o->has_arg != copt_no_argument) { - if (argv[opt->_i][j] == '=') - opt->arg = &argv[opt->_i][j + 1]; - else if (argv[opt->_i][j] == '\0' && opt->_i < argc - 1 && (o->has_arg == copt_required_argument || - argv[opt->_i + 1][0] != '-')) - opt->arg = argv[++opt->_i]; - else if (o->has_arg == copt_required_argument) - optc = ':'; /* missing option argument */ - } - } - } - } else { /* a short option */ - const char *p; - if (opt->_pos == 0) opt->_pos = 1; - optc = opt->opt = argv[opt->_i][opt->_pos++]; - opt->_faulty[1] = optc, opt->faulty = opt->_faulty; - p = strchr((char *) shortopts, optc); - if (p == 0) { - optc = '?'; /* unknown option */ - } else if (p[1] == ':') { - if (argv[opt->_i][opt->_pos] != '\0') - opt->arg = &argv[opt->_i][opt->_pos]; - else if (opt->_i < argc - 1 && (p[2] != ':' || argv[opt->_i + 1][0] != '-')) - opt->arg = argv[++opt->_i]; - else if (p[2] != ':') - optc = ':'; - opt->_pos = -1; - } - } - if (opt->_pos < 0 || argv[opt->_i][opt->_pos] == 0) { - ++opt->_i, opt->_pos = 0; - if (opt->_nargs > 0) /* permute */ - for (j = i0; j < opt->_i; ++j) - _copt_permute(argv, j, opt->_nargs); - } - opt->ind = opt->_i - opt->_nargs; - return optc; -} - -#endif diff --git a/stc/crandom.h b/stc/crandom.h index af7070bc..bc257f35 100644 --- a/stc/crandom.h +++ b/stc/crandom.h @@ -27,192 +27,82 @@ #include "cdefs.h" #include -/* - * PCG32 random number generator: https://www.pcg-random.org/index.html - */ -typedef struct {uint64_t state; uint64_t inc;} crandom32_t; +/* ------------- CRand32: PCG32 -------------- */ -STC_INLINE uint32_t crandom32(crandom32_t* rng) -{ - uint64_t old = rng->state; - rng->state = old * 6364136223846793005ull + rng->inc; - uint32_t xos = ((old >> 18u) ^ old) >> 27u; - uint32_t rot = old >> 59u; - return (xos >> rot) | (xos << ((-rot) & 31)); -} +typedef struct {uint64_t state, inc;} CRand32; -/* float random int number in range [0, 1). Note: 23 bit resolution. */ -STC_INLINE float crandom32f(crandom32_t* rng) { - union {uint32_t i; float f;} u = {0x3F800000u | (crandom32(rng) >> 9)}; - return u.f - 1.0f; -} +STC_API CRand32 crand32_init2(uint64_t seed, uint64_t seq); -/* Uniform random number in range [0, bound) */ -STC_INLINE uint32_t crandom32b(crandom32_t* rng, uint32_t bound) { - return (uint32_t) (((uint64_t) crandom32(rng) * bound) >> 32); +STC_INLINE CRand32 crand32_init(uint64_t seed) { + return crand32_init2(seed, seed); } +STC_API uint32_t crand32_gen(CRand32* rng); -STC_INLINE crandom32_t crandom32_init2(uint64_t seed, uint64_t seq) { - crandom32_t rng = {0u, (seq << 1u) | 1u}; /* inc must be odd */ - crandom32(&rng); - rng.state += seed; - crandom32(&rng); - return rng; +/* Uniform random number in range [0, bound) */ +STC_INLINE uint32_t crand32_genBounded(CRand32* rng, uint32_t bound) { + return (uint32_t) (((uint64_t) crand32_gen(rng) * bound) >> 32); } -STC_INLINE crandom32_t crandom32_init(uint64_t seed) { - return crandom32_init2(seed, seed); +/* float random int number in range [0, 1). Note: 23 bit resolution. */ +STC_INLINE float crand32_genReal(CRand32* rng) { + union {uint32_t i; float f;} u = {0x3F800000u | (crand32_gen(rng) >> 9)}; + return u.f - 1.0f; } +/* ------------- CRand64: SFC64 -------------- */ -/* Rotate bits left */ -STC_INLINE uint64_t c_rotateLeft64(uint64_t x, int bits) { - return (x << bits) | (x >> (64 - bits)); -} +typedef struct {uint64_t state[4];} CRand64; -/* - * SFC64 random number generator: http://pracrand.sourceforge.net - */ -typedef struct {uint64_t state[4];} crandom64_t; +STC_API CRand64 crand64_init(const uint64_t seed); -STC_API uint64_t crandom64(crandom64_t* rng) { - enum {LR=24, RS=11, LS=3}; - uint64_t *s = rng->state; - const uint64_t result = s[0] + s[1] + s[3]++; - s[0] = s[1] ^ (s[1] >> RS); - s[1] = s[2] + (s[2] << LS); - s[2] = c_rotateLeft64(s[2], LR) + result; - return result; -} +STC_API uint64_t crand64_gen(CRand64* rng); /* float64 random int number in range [0, 1), 52 bit resolution */ -STC_INLINE double crandom64f(crandom64_t* rng) { - union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (crandom64(rng) >> 12)}; +STC_INLINE double crand64_genReal(CRand64* rng) { + union {uint64_t i; double f;} u = {0x3FF0000000000000ull | (crand64_gen(rng) >> 12)}; return u.f - 1.0; } -STC_API crandom64_t crandom64_init(const uint64_t seed) { - crandom64_t state = {{seed, seed, seed, 1}}; - for (int i = 0; i < 12; ++i) crandom64(&state); - return state; -} -/* - * SipHash implementation. - */ -#if defined(_WIN32) || (defined(__BYTE_ORDER__) && __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__) - STC_INLINE uint64_t c_le64ToHost(uint64_t x) { return x; } -#elif defined(__APPLE__) - #include - STC_INLINE uint64_t c_le64ToHost(uint64_t x) { return OSSwapLittleToHostInt64(x); } -#elif defined(__FreeBSD__) || defined(__NetBSD__) || defined(__OpenBSD__) || defined(__DragonFly__) - #include - STC_INLINE uint64_t c_le64ToHost(uint64_t x) { return letoh64(x); } -#elif defined(__linux__) || defined(__CYGWIN__) || defined(__GNUC__) || defined(__GNU_LIBRARY__) - #include - STC_INLINE uint64_t c_le64ToHost(uint64_t x) { return le64toh(x); } -#endif +#if !defined(STC_HEADER) || defined(STC_IMPLEMENTATION) -typedef struct siphash_t { - uint64_t v[4], padding; - size_t length; - int c, d; -} siphash_t; - -/* c=2, d=4 or c=1, d=3 */ -STC_INLINE void siphash_init_c_d(siphash_t* s, const uint64_t key[2], const int c, const int d) { - s->c = c; - s->d = d; - s->length = 0; - s->padding = 0; - s->v[0] = key[0] ^ 0x736f6d6570736575; - s->v[1] = key[1] ^ 0x646f72616e646f6d; - s->v[2] = key[0] ^ 0x6c7967656e657261; - s->v[3] = key[1] ^ 0x7465646279746573; -} +/* PCG32 random number generator: https://www.pcg-random.org/index.html */ -/* default init 2-4 */ -STC_API siphash_t siphash_init(const uint64_t key[2]) { - siphash_t state; - siphash_init_c_d(&state, key, 2, 4); - return state; +STC_API CRand32 crand32_init2(uint64_t seed, uint64_t seq) { + CRand32 rng = {0u, (seq << 1u) | 1u}; /* inc must be odd */ + crand32_gen(&rng); + rng.state += seed; + crand32_gen(&rng); + return rng; } -#define _siphash_halfRound(i, j, a, b, c, d) \ - (a += b, \ - c += d, \ - b = c_rotateLeft64(b, i) ^ a, \ - d = c_rotateLeft64(d, j) ^ c, \ - a = c_rotateLeft64(a, 32)) - -#define _siphash_compress(rounds, v) \ - for (int r = 0; r < rounds; ++r) { \ - _siphash_halfRound(13, 16, v[0], v[1], v[2], v[3]); \ - _siphash_halfRound(17, 21, v[2], v[1], v[0], v[3]); \ - } - -#define _siphash_digest(rounds, v, m) { \ - const uint64_t _m = m; \ - v[3] ^= _m; \ - _siphash_compress(rounds, v); \ - v[0] ^= _m; \ - } - -STC_API void siphash_update(siphash_t* s, const void* bytes, size_t size) { - union { const uint8_t* u8; const uint64_t* u64; } in; - in.u8 = (const uint8_t*) bytes; - size_t offset = s->length & 7; - uint64_t *v = s->v; - s->length += size; - - if (offset) { - const size_t end = offset + size; - size -= 8 - offset; - while (offset < end && offset < 8) { - s->padding |= ((uint64_t) *in.u8++) << (offset++ << 3); - } - if (end < 8) return; - - _siphash_digest(s->c, v, s->padding); - s->padding = 0; - } - size_t n_words = size >> 3; - uint64_t m; - - while (n_words--) { - memcpy(&m, in.u64++, 8); - _siphash_digest(s->c, v, c_le64ToHost(m)); - } - switch (s->length & 7) { - case 7: s->padding |= ((uint64_t) in.u8[6]) << 48; - case 6: s->padding |= ((uint64_t) in.u8[5]) << 40; - case 5: s->padding |= ((uint64_t) in.u8[4]) << 32; - case 4: s->padding |= ((uint64_t) in.u8[3]) << 24; - case 3: s->padding |= ((uint64_t) in.u8[2]) << 16; - case 2: s->padding |= ((uint64_t) in.u8[1]) << 8; - case 1: s->padding |= ((uint64_t) in.u8[0]); - } +STC_API uint32_t crand32_gen(CRand32* rng) { + uint64_t old = rng->state; + rng->state = old * 6364136223846793005ull + rng->inc; + uint32_t xos = ((old >> 18u) ^ old) >> 27u; + uint32_t rot = old >> 59u; + return (xos >> rot) | (xos << ((-rot) & 31)); } -STC_API uint64_t siphash_finalize(siphash_t* s) { - uint64_t *v = s->v; - _siphash_digest(s->c, v, s->padding | (s->length << 56)); - v[2] ^= 0xff; - _siphash_compress(s->d, v); - return v[0] ^ v[1] ^ v[2] ^ v[3]; -} -/* c=2, d=4 or c=1, d=3 */ -STC_API uint64_t siphash_hash_c_d(const uint64_t key[2], const void* bytes, const uint64_t size, const int c, const int d) { - siphash_t state; - siphash_init_c_d(&state, key, c, d); - siphash_update(&state, bytes, size); - return siphash_finalize(&state); +/* SFC64 random number generator: http://pracrand.sourceforge.net */ + +STC_API CRand64 crand64_init(const uint64_t seed) { + CRand64 state = {{seed, seed, seed, 1}}; + for (int i = 0; i < 12; ++i) crand64_gen(&state); + return state; } -/* default hash 2-4 */ -STC_INLINE uint64_t siphash_hash(const uint64_t key[2], const void* bytes, const uint64_t size) { - return siphash_hash_c_d(key, bytes, size, 2, 4); +STC_API uint64_t crand64_gen(CRand64* rng) { + enum {LR=24, RS=11, LS=3}; + uint64_t *s = rng->state; + const uint64_t result = s[0] + s[1] + s[3]++; + s[0] = s[1] ^ (s[1] >> RS); + s[1] = s[2] + (s[2] << LS); + s[2] = (s[2] << LR) | (s[2] >> (64 - LR)) + result; + return result; } #endif + +#endif -- cgit v1.2.3