From 29365e9144e7e5633ed4d9b83ed7eafa71768a02 Mon Sep 17 00:00:00 2001 From: Reini Urban Date: Sun, 22 Nov 2015 16:09:04 +0100 Subject: [PATCH] replace 32bit mersenne twister with pcg and added bounds argument for rand and num rand. See GH #99 But it is 6.8% slower than mt: 3m11.529s vs 2m58.521s for binarytrees.pn --- ChangeLog | 2 + Makefile | 4 +- core/mt19937ar.c | 184 ----------------------------------------- core/number.c | 4 +- core/objmodel.c | 2 +- core/pcg.c | 211 +++++++++++++++++++++++++++++++++++++++++++++++ core/potion.h | 2 +- 7 files changed, 219 insertions(+), 190 deletions(-) delete mode 100644 core/mt19937ar.c create mode 100644 core/pcg.c diff --git a/ChangeLog b/ChangeLog index 834ee44c..edaf36d9 100644 --- a/ChangeLog +++ b/ChangeLog @@ -5,6 +5,8 @@ v0.4 * Fix installation from source without git * Fixed various SEGVs and minor errors and security issues. * Fixed cmp for string (Peter Arthur) + * Replaced 32bit mersenne twister with pcg. + * Added bounds argument for rand and num rand. v0.3 Sat Oct 17 10:42:49 CEST 2015 rurban v0.3-g7963cbc-1264 diff --git a/Makefile b/Makefile index 88f94e82..15dd4495 100644 --- a/Makefile +++ b/Makefile @@ -4,7 +4,7 @@ examples release deploy install grammar doxygen website testable .NOTPARALLEL: test -SRC = core/asm.c core/ast.c core/compile.c core/contrib.c core/gc.c core/internal.c core/lick.c core/mt19937ar.c core/number.c core/objmodel.c core/primitive.c core/string.c core/syntax.c core/table.c core/vm.c +SRC = core/asm.c core/ast.c core/compile.c core/contrib.c core/gc.c core/internal.c core/lick.c core/pcg.c core/number.c core/objmodel.c core/primitive.c core/string.c core/syntax.c core/table.c core/vm.c PLIBS = readline buffile aio database PLIBS_SRC = lib/aio.c lib/buffile.c lib/database.c lib/readline/readline.c lib/readline/linenoise.c GREGCFLAGS = -O3 -DNDEBUG @@ -202,7 +202,7 @@ $(foreach o,${OBJS},core/internal.${o} ): core/internal.c core/potion.h core/con $(foreach o,${OBJS},core/lick.${o} ): core/lick.c core/potion.h core/config.h core/internal.h $(foreach o,${OBJS},core/load.${o} ): core/load.c core/potion.h core/config.h core/internal.h \ core/table.h core/khash.h -$(foreach o,${OBJS},core/mt19937ar.${o} ): core/mt19937ar.c core/potion.h core/config.h +$(foreach o,${OBJS},core/pcg.${o} ): core/pcg.c core/potion.h core/config.h $(foreach o,${OBJS},core/number.${o} ): core/number.c core/potion.h core/config.h core/internal.h $(foreach o,${OBJS},core/objmodel.${o} ): core/objmodel.c core/potion.h core/config.h core/internal.h \ core/khash.h core/table.h core/asm.h diff --git a/core/mt19937ar.c b/core/mt19937ar.c deleted file mode 100644 index 516fa787..00000000 --- a/core/mt19937ar.c +++ /dev/null @@ -1,184 +0,0 @@ -/**\file mt19937ar.c - random numbers (mersenne twister). As Lobby (global long) or PNInteger (0-1.0) - - A C-program for MT19937, with initialization improved 2002/2/10. - Coded by Takuji Nishimura and Makoto Matsumoto. - This is a faster version by taking Shawn Cokus's optimization, - Matthe Bellew's simplification, Isaku Wada's real version. - - Before using, initialize the state by using init_genrand(seed) - or init_by_array(init_key, key_length). - - Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, - All rights reserved. - - Redistribution and use in source and binary forms, with or without - modification, are permitted provided that the following conditions - are met: - - 1. Redistributions of source code must retain the above copyright - notice, this list of conditions and the following disclaimer. - - 2. Redistributions in binary form must reproduce the above copyright - notice, this list of conditions and the following disclaimer in the - documentation and/or other materials provided with the distribution. - - 3. The names of its contributors may not be used to endorse or promote - products derived from this software without specific prior written - permission. - - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS - "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT - LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR - A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR - CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR - PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF - LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING - NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS - SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - - - Any feedback is very welcome. - http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html - email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) -*/ - -#include -#include "potion.h" - -/* Period parameters */ -#define N 624 -#define M 397 -#define MATRIX_A 0x9908b0dfUL /* constant vector a */ -#define UMASK 0x80000000UL /* most significant w-r bits */ -#define LMASK 0x7fffffffUL /* least significant r bits */ -#define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) ) -#define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL)) - -static unsigned long state[N]; /* the array for the state vector */ -static int left = 1; -static int initf = 0; -static unsigned long *next; - -/* initializes state[N] with a seed */ -void init_genrand(unsigned long s) { - int j; - state[0]= s & 0xffffffffUL; - for (j=1; j> 30)) + j); - /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ - /* In the previous versions, MSBs of the seed affect */ - /* only MSBs of the array state[]. */ - /* 2002/01/09 modified by Makoto Matsumoto */ - state[j] &= 0xffffffffUL; /* for >32 bit machines */ - } - left = 1; initf = 1; -} - -/* initialize by an array with array-length */ -/* init_key is the array for initializing keys */ -/* key_length is its length */ -/* slight change for C++, 2004/2/26 */ -void init_by_array(unsigned long init_key[], int key_length) { - int i, j, k; - init_genrand(19650218UL); - i=1; j=0; - k = (N>key_length ? N : key_length); - for (; k; k--) { - state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL)) - + init_key[j] + j; /* non linear */ - state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ - i++; j++; - if (i>=N) { state[0] = state[N-1]; i=1; } - if (j>=key_length) j=0; - } - for (k=N-1; k; k--) { - state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL)) - - i; /* non linear */ - state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ - i++; - if (i>=N) { state[0] = state[N-1]; i=1; } - } - - state[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ - left = 1; initf = 1; -} - -static void next_state(void) { - unsigned long *p=state; - int j; - - /* if init_genrand() has not been called, */ - /* a default initial seed is used */ - if (initf==0) init_genrand(5489UL); - - left = N; - next = state; - - for (j=N-M+1; --j; p++) { - /*@ assert Value: mem_access: \valid_read(p+1); */ - /*@ assert Value: mem_access: \valid_read(p+397); */ - *p = p[M] ^ TWIST(p[0], p[1]); - } - - for (j=M; --j; p++) { - /*@ assert Value: mem_access: \valid_read(p+1); */ - /*@ assert Value: mem_access: \valid_read(p+(int)(397-624)); */ - *p = p[M-N] ^ TWIST(p[0], p[1]); - } - - /*@ assert Value: mem_access: \valid(p); */ - /*@ assert Value: mem_access: \valid_read(p+(int)(397-624)); */ - *p = p[M-N] ^ TWIST(p[0], state[0]); -} - -/** generates a random number on [0,0xffffffff]-interval */ -unsigned long potion_rand_int(void) { - unsigned long y; - - if (--left == 0) next_state(); - y = *next++; - - /* Tempering */ - y ^= (y >> 11); - y ^= (y << 7) & 0x9d2c5680UL; - y ^= (y << 15) & 0xefc60000UL; - y ^= (y >> 18); - - return y; -} - -/** generates a random number on [0,1) with 53-bit resolution*/ -double potion_rand_double(void) { - unsigned long a=potion_rand_int()>>5, b=potion_rand_int()>>6; - return(a*67108864.0+b)*(1.0/9007199254740992.0); -} -/**\memberof Lobby - "srand" initialize random seed - \param seed PNInteger - \return Lobby (usually unused) - \sa potion_rand. */ -PN potion_srand(Potion *P, PN cl, PN self, PN seed) { - init_genrand(PN_INT(seed)); - return self; -} - -/**\memberof Lobby - "rand" generate random ulong number - \code rand #=> xxxxxx \endcode - \return PNDouble in [0,0xffffffff]-interval - \sa potion_num_rand for double, potion_srand. */ -PN potion_rand(Potion *P, PN cl, PN self) { - return PN_NUM(potion_rand_int()); -} - -/**\memberof PNInteger - "rand" generate random float number - \code 1 rand #=> 0.xxxxxx \endcode - \return PNDouble in [0.0,0.1]-interval - \sa potion_rand for long, potion_srand. */ -PN potion_num_rand(Potion *P, PN cl, PN self) { - return PN_NUM(potion_rand_double()); -} diff --git a/core/number.c b/core/number.c index 2d674f97..68fdbcb5 100644 --- a/core/number.c +++ b/core/number.c @@ -469,7 +469,7 @@ void potion_num_init(Potion *P) { potion_method(dbl_vt, "/", potion_dbl_div, "value=D"); potion_method(dbl_vt, "abs", potion_dbl_abs, 0); potion_method(dbl_vt, "cmp", potion_dbl_cmp, "value=D"); - //potion_method(dbl_vt, "rand", potion_dbl_rand, 0); + //potion_method(dbl_vt, "rand", potion_dbl_rand, "|bound=D"); // optimized integer-only methods, for both operands potion_method(int_vt, "+", potion_int_add, "value=I"); potion_method(int_vt, "-", potion_int_sub, "value=I"); @@ -482,7 +482,7 @@ void potion_num_init(Potion *P) { potion_method(int_vt, "chr", potion_int_chr, 0); potion_method(int_vt, "abs", potion_int_abs, 0); potion_method(int_vt, "cmp", potion_int_cmp, "value=I"); - //potion_method(int_vt, "rand", potion_int_rand, 0); + //potion_method(int_vt, "rand", potion_int_rand, "|bound=I"); potion_method(num_vt, "step", potion_int_step, "end=N,step=N,block=&"); potion_method(num_vt, "times", potion_int_times, "block=&"); potion_method(num_vt, "to", potion_int_to, "end=N,block=&"); diff --git a/core/objmodel.c b/core/objmodel.c index a25e7cdb..d36d6044 100644 --- a/core/objmodel.c +++ b/core/objmodel.c @@ -748,7 +748,7 @@ void potion_lobby_init(Potion *P) { potion_method(P->lobby, "kind", potion_lobby_kind, 0); potion_method(P->lobby, "isa?", potion_lobby_isa, "value=o"); potion_method(P->lobby, "srand", potion_srand, "seed=N"); - potion_method(P->lobby, "rand", potion_rand, 0); + potion_method(P->lobby, "rand", potion_rand, "|bound=N"); potion_method(P->lobby, "self", potion_lobby_self, 0); potion_method(P->lobby, "string", potion_lobby_string, 0); potion_method(P->lobby, "can", potion_lobby_can, "method=S"); diff --git a/core/pcg.c b/core/pcg.c new file mode 100644 index 00000000..dd0f8922 --- /dev/null +++ b/core/pcg.c @@ -0,0 +1,211 @@ +/* + * PCG Random Number Generation for C. Based on pcg_basic.c + * + * Copyright 2014 Melissa O'Neill + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * For additional information about the PCG random number generation scheme, + * including its license and other licensing options, visit + * + * http://www.pcg-random.org + */ + +/* + * This code is derived from the full C implementation, which is in turn + * derived from the canonical C++ PCG implementation. The C++ version + * has many additional features and is preferable if you can use C++ in + * your project. + */ + +#include "potion.h" +#include +#ifndef _WIN32 +#include +#else +#define WIN32_LEAN_AND_MEAN +#define NOGDI +#include /* INT_MAX */ +typedef unsigned __int32 uint32_t; +typedef unsigned __int64 uint64_t; +#define inline __inline +#endif + + +struct pcg_state_setseq_64 { /* Internals are *Private*. */ + uint64_t state; /* RNG state. All values are possible. */ + uint64_t inc; /* Controls which RNG sequence (stream) is */ + /* selected. Must *always* be odd. */ +}; +typedef struct pcg_state_setseq_64 pcg32_random_t; + +/* If you *must* statically initialize it, here's one. */ + +#define PCG32_INITIALIZER { 0x853c49e6748fea9bULL, 0xda3e39cb94b95bdbULL } + +/* state for global RNGs */ + +static pcg32_random_t pcg32_global = PCG32_INITIALIZER; + +/* pcg32_srandom(initstate, initseq) + * Seed the rng. Specified in two parts, state initializer and a + * sequence selection constant (a.k.a. stream id) */ + +void pcg32_srandom(uint64_t initstate, uint64_t initseq); + +/* pcg32_random() + * Generate a uniformly distributed 32-bit random number */ + +uint32_t pcg32_random(void); +static uint32_t pcg32_random_r(pcg32_random_t* rng); + +/* pcg32_boundedrand(bound) + * Generate a uniformly distributed number, r, where 0 <= r < bound */ + +uint32_t pcg32_boundedrand(uint32_t bound); + +/* pcg32_srandom(initstate, initseq) + * Seed the rng. Specified in two parts, state initializer and a + * sequence selection constant (a.k.a. stream id) */ + +static void pcg32_srandom_r(pcg32_random_t* rng, uint64_t initstate, uint64_t initseq) +{ + rng->state = 0U; + rng->inc = (initseq << 1u) | 1u; + pcg32_random_r(rng); + rng->state += initstate; + pcg32_random_r(rng); +} + +void pcg32_srandom(uint64_t seed, uint64_t seq) +{ + pcg32_srandom_r(&pcg32_global, seed, seq); +} + +/* pcg32_random() + * pcg32_random_r(rng) + * Generate a uniformly distributed 32-bit random number */ + +static uint32_t pcg32_random_r(pcg32_random_t* rng) +{ + uint64_t oldstate = rng->state; + rng->state = oldstate * 6364136223846793005ULL + rng->inc; + uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u; + uint32_t rot = oldstate >> 59u; + return (xorshifted >> rot) | (xorshifted << ((-rot) & 31)); +} + +uint32_t pcg32_random(void) +{ + return pcg32_random_r(&pcg32_global); +} + + +/* pcg32_boundedrand(bound): + * pcg32_boundedrand_r(rng, bound): + * Generate a uniformly distributed number, r, where 0 <= r < bound */ + +static uint32_t pcg32_boundedrand_r(pcg32_random_t* rng, uint32_t bound) +{ + /* To avoid bias, we need to make the range of the RNG a multiple of + * bound, which we do by dropping output less than a threshold. + * A naive scheme to calculate the threshold would be to do + * + * uint32_t threshold = 0x100000000ull % bound; + * + * but 64-bit div/mod is slower than 32-bit div/mod (especially on + * 32-bit platforms). In essence, we do + * + * uint32_t threshold = (0x100000000ull-bound) % bound; + * + * because this version will calculate the same modulus, but the LHS + * value is less than 2^32. */ + + uint32_t threshold = -bound % bound; + + /* Uniformity guarantees that this loop will terminate. In practice, it + * should usually terminate quickly; on average (assuming all bounds are + * equally likely), 82.25% of the time, we can expect it to require just + * one iteration. In the worst case, someone passes a bound of 2^31 + 1 + * (i.e., 2147483649), which invalidates almost 50% of the range. In + * practice, bounds are typically small and only a tiny amount of the range + * is eliminated. */ + for (;;) { + uint32_t r = pcg32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +uint32_t pcg32_boundedrand(uint32_t bound) +{ + return pcg32_boundedrand_r(&pcg32_global, bound); +} + +/*********************************************************/ + +/** generates a random number on [0,0xffffffff]-interval */ +unsigned long potion_rand_int(void) { + return pcg32_random_r(&pcg32_global); +} + +/** generates a random number on [0, 1.0) with 53-bit resolution*/ +double potion_rand_double(void) { + unsigned long a = potion_rand_int()>>5, b = potion_rand_int()>>6; + return(a*67108864.0+b)*(1.0/9007199254740992.0); +} + +/**\memberof Lobby + "srand" initialize random seed + \param seed PNInteger + \return Lobby (usually unused) + \sa potion_rand. */ +PN potion_srand(Potion *P, PN cl, PN self, PN seed) { + pcg32_srandom_r(&pcg32_global, PN_INT(seed), 1); // TODO: a stream id + return self; +} + +/**\memberof Lobby + "rand" generate random ulong number + \code rand #=> xxxxxx \endcode + \param bound PNInteger, default: 0xffffffff + \return PNInteger in [0, bound) interval + \sa potion_num_rand for double, potion_srand. */ +PN potion_rand(Potion *P, PN cl, PN self, PN bound) { + return bound ? PN_NUM(pcg32_boundedrand(PN_INT(bound))) + : PN_NUM(potion_rand_int()); +} + +/**\memberof PNNumber + "rand" generate random float number + \code 0 rand #=> 0.xxxxxx \endcode + \code 1 rand #=> 0.xxxxxx \endcode + \code 2.0 rand #=> [0-1].xxxxxx \endcode + \return PNDouble in [0.0, self) interval + \sa potion_rand for long, potion_srand. */ +PN potion_num_rand(Potion *P, PN cl, PN self) { + double bound = PN_DBL(self); + if (bound == 0.0 || bound == 1.0) { + return potion_double(P, potion_rand_double()); + } else { + if (bound == trunc(bound)) // 2.0 => [0-1] + dblrand + return potion_double(P, pcg32_boundedrand((uint32_t)bound-1) + potion_rand_double()); + else { + double x; + do { + x = pcg32_boundedrand((uint32_t)bound) + potion_rand_double(); + } while (x >= bound); + return potion_double(P, x); + } + } +} diff --git a/core/potion.h b/core/potion.h index 008754eb..0d0deaab 100644 --- a/core/potion.h +++ b/core/potion.h @@ -859,7 +859,7 @@ PN potion_strtod(Potion *, char *, int); PN potion_num_pow(Potion *, PN, PN, PN); PN potion_num_rand(Potion *, PN, PN); PN potion_srand(Potion *, PN, PN, PN); -PN potion_rand(Potion *, PN, PN); +PN potion_rand(Potion *, PN, PN, PN); PN potion_sig_at(Potion *, PN, int); PN potion_sig_name_at(Potion *, PN, int); int potion_sig_arity(Potion *, PN);