From acd71d910719815f43caeac7a08be8f1dae7a609 Mon Sep 17 00:00:00 2001 From: Chris Jefferson Date: Fri, 28 Sep 2018 11:10:34 +0100 Subject: [PATCH] Add PCG32 random number generator This random generator produces very high quality random numbers. It also has the major advantage that it's initialisation and state are both very small and simple, meaning it is very low cost to create many instances of this generator --- lib/random.gd | 14 +++- lib/random.gi | 53 ++++++++++++++- src/intfuncs.c | 163 ++++++++++++++++++++++++++++++++++++++++++++++- tst/testrandom.g | 6 +- 4 files changed, 230 insertions(+), 6 deletions(-) diff --git a/lib/random.gd b/lib/random.gd index 9e5a9b2787..3103a857d4 100644 --- a/lib/random.gd +++ b/lib/random.gd @@ -146,11 +146,12 @@ DeclareOperation( "Init", [IsRandomSource, IsObject] ); ## ## ## +## ## ## ## ## -## Currently, the ⪆ library provides three types of random sources, +## Currently, the ⪆ library provides four types of random sources, ## distinguished by the three listed categories. ##

## are random sources which use a fast @@ -160,8 +161,14 @@ DeclareOperation( "Init", [IsRandomSource, IsObject] ); ## and the origin of the code used in the ⪆ kernel, see: ## http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html. ##

-## Use the Mersenne twister if possible, in particular for generating many -## large random integers. +## The Mersenne twister is a good standard choice, in particular for +## generating many large random integers. +##

+## are random sources which use the PCG +## random generator. These generators achieve excellent statistical quality +## and also have the advantage that their internal state is very small +## (two 64-bit integers), so can be created very quickly. PCG32 is documented +## at http://www.pcg-random.org/. ##

## There is also a predefined global random source ## which is used by most of the library @@ -187,6 +194,7 @@ DeclareOperation( "Init", [IsRandomSource, IsObject] ); DeclareCategory("IsGlobalRandomSource", IsRandomSource); DeclareCategory("IsGAPRandomSource", IsRandomSource); DeclareCategory("IsMersenneTwister", IsRandomSource); +DeclareCategory("IsPCG32RandomSource", IsRandomSource); if IsHPCGAP then MakeThreadLocal( "GlobalRandomSource" ); diff --git a/lib/random.gi b/lib/random.gi index db56ff281d..4f46118ad5 100644 --- a/lib/random.gi +++ b/lib/random.gi @@ -31,7 +31,7 @@ end); # Generic fallback methods, such that it is sufficient to install Random for # lists or for pairs of integers. -InstallMethod(Random, [IsRandomSource, IsInt, IsInt], function(rs, a, b) +BindGlobal("FALLBACK_GENERATE_RANDOM_INT", function(rs, a, b) local d, x, r, y; d := b - a; if d < 0 then @@ -54,6 +54,8 @@ InstallMethod(Random, [IsRandomSource, IsInt, IsInt], function(rs, a, b) fi; end); +InstallMethod(Random, [IsRandomSource, IsInt, IsInt], FALLBACK_GENERATE_RANDOM_INT); + InstallMethod(Random, [IsRandomSource, IsList], function(rs, list) return list[Random(rs, 1, Length(list))]; end); @@ -159,6 +161,55 @@ InstallMethod(Random, [IsGAPRandomSource, IsList], function(rs, list) fi; end); +############################################################################ +## PCG32 Random Source +## + +InstallMethod(Init, [IsPCG32RandomSource, IsObject], function(rs, seed) + local initstate, initseq; + initstate := 1; + initseq := 0; + if IsInt(seed) then + rs!.state := PCG32_INIT_FROM_SEED(seed, 0); + else + rs!.state := PCG32_INIT_FROM_STATE(seed[1], seed[2]); + fi; + return rs; +end); + +InstallMethod(State, [IsPCG32RandomSource], function(rs) + return PCG32_GET_STATE(rs!.state); +end); + +InstallMethod(Reset, [IsPCG32RandomSource, IsObject], function(rs, seed) + local old; + old := State(rs); + Init(rs, seed); + return old; +end); + +InstallMethod(Random, [IsPCG32RandomSource, IsList], function(rs, list) + local val; + if Length(list) < 2^32-1 then + return list[PCG32_BOUNDEDRAND_R(rs!.state, Length(list))+1]; + else + return list[Random(rs, 1, Length(list))]; + fi; +end); + +InstallMethod(Random, [IsPCG32RandomSource, IsInt, IsInt], function(rs,a,b) + local d, x, r, y; + d := b - a; + if d < 0 then + return fail; + elif a = b then + return a; + elif d < 2^32-1 then + return a + PCG32_BOUNDEDRAND_R(rs!.state, d+1); + else + return FALLBACK_GENERATE_RANDOM_INT(rs,a,b); + fi; +end); ############################################################################## ## Random source using the Mersenne twister kernel functions. diff --git a/src/intfuncs.c b/src/intfuncs.c index d7a45eea77..8c79a013ed 100644 --- a/src/intfuncs.c +++ b/src/intfuncs.c @@ -18,6 +18,7 @@ #include "bool.h" #include "calls.h" #include "error.h" +#include "integer.h" #include "lists.h" #include "modules.h" #include "plist.h" @@ -26,6 +27,159 @@ #include "stringobj.h" +/**************************************************************************** +** +** * * * * * * * "PCG" random numbers * * * * * * * * * * * * * +** +** This next function implements the stepping function of the 'PCG' +** random number library. This library produces high quality random numbers +** which pass many statistical tests. They also have a very small state which +** can be quickly initalised +** +** The following three functions are under the following original license +** *Really* minimal PCG32 code / (c) 2014 M.E. O'Neill / pcg-random.org +** Licensed under Apache License 2.0 (NO WARRANTY, etc. see website) +*/ + +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; + +// 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)); +} + +// pcg32_srandom(initstate, initseq) +// pcg32_srandom_r(rng, 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); +} + +// 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; + } +} + +// End of external PCG32 code + + +static Obj TYPE_KERNEL_OBJECT; + +Obj FuncPCG32_INIT_FROM_SEED(Obj self, Obj initstateObj, Obj initseqObj) +{ + uint64_t initstate = UInt8_ObjInt(initstateObj); + uint64_t initseq = UInt8_ObjInt(initseqObj); + + Obj state = NewBag(T_DATOBJ, sizeof(Obj) + sizeof(pcg32_random_t)); + SET_TYPE_DATOBJ(state, TYPE_KERNEL_OBJECT); + pcg32_srandom_r((pcg32_random_t *)(1 + ADDR_OBJ(state)), initstate, + initseq); + return state; +} + +// This function assumes we have an existing state we are reinitalising from. +// It will reject if initseq is not odd. +Obj FuncPCG32_INIT_FROM_STATE(Obj self, Obj initstateObj, Obj initseqObj) +{ + uint64_t initstate = UInt8_ObjInt(initstateObj); + uint64_t initseq = UInt8_ObjInt(initseqObj); + + if (initseq % 2 == 0) { + ErrorMayQuit("Illegal PCG32 state: 2nd argument must be odd", 0, 0); + } + + Obj state = NewBag(T_DATOBJ, sizeof(Obj) + sizeof(pcg32_random_t)); + SET_TYPE_DATOBJ(state, TYPE_KERNEL_OBJECT); + pcg32_random_t * pcg32_state = (pcg32_random_t *)(1 + ADDR_OBJ(state)); + pcg32_state->state = initstate; + pcg32_state->inc = initseq; + return state; +} + +Obj FuncPCG32_RANDOM_R(Obj self, Obj state) +{ + uint32_t val = pcg32_random_r((pcg32_random_t *)(1 + ADDR_OBJ(state))); + return ObjInt_UInt(val); +} + +Obj FuncPCG32_BOUNDEDRAND_R(Obj self, Obj state, Obj boundObj) +{ + UInt8 bound = UInt8_ObjInt(boundObj); + if (bound > UINT32_MAX) { + ErrorMayQuit("Invalid bound in PCG32_BOUNBDEDRAND_R", 0, 0); + } + + uint32_t val = + pcg32_boundedrand_r((pcg32_random_t *)(1 + ADDR_OBJ(state)), bound); + return ObjInt_UInt(val); +} + +Obj FuncPCG32_GET_STATE(Obj self, Obj stateObj) +{ + pcg32_random_t * pcg32_state = (pcg32_random_t *)(1 + ADDR_OBJ(stateObj)); + Obj state = ObjInt_UInt8(pcg32_state->state); + Obj inc = ObjInt_UInt8(pcg32_state->inc); + Obj list = NEW_PLIST(T_PLIST, 2); + SET_LEN_PLIST(list, 2); + SET_ELM_PLIST(list, 1, state); + SET_ELM_PLIST(list, 2, inc); + CHANGED_BAG(list); + return list; +} + /**************************************************************************** ** ** * * * * * * * "Mersenne twister" random numbers * * * * * * * * * * * * * @@ -730,13 +884,18 @@ Obj FuncMAKE_BITFIELDS(Obj self, Obj widths) ** *V GVarFuncs . . . . . . . . . . . . . . . . . . list of functions to export */ -static StructGVarFunc GVarFuncs [] = { +static StructGVarFunc GVarFuncs[] = { GVAR_FUNC(HASHKEY_BAG, 4, "obj, int,int,int"), GVAR_FUNC(InitRandomMT, 1, "initstr"), GVAR_FUNC(MAKE_BITFIELDS, -1, "widths"), GVAR_FUNC(BUILD_BITFIELDS, -2, "widths, vals"), + GVAR_FUNC(PCG32_INIT_FROM_SEED, 2, "initstate, initseq"), + GVAR_FUNC(PCG32_INIT_FROM_STATE, 2, "initstate, initseq"), + GVAR_FUNC(PCG32_RANDOM_R, 1, "state"), + GVAR_FUNC(PCG32_BOUNDEDRAND_R, 2, "state, bound"), + GVAR_FUNC(PCG32_GET_STATE, 1, "state"), { 0, 0, 0, 0, 0 } }; @@ -758,6 +917,8 @@ static Int InitKernel ( /* init filters and functions */ InitHdlrFuncsFromTable( GVarFuncs ); + ImportGVarFromLibrary("TYPE_KERNEL_OBJECT", &TYPE_KERNEL_OBJECT); + return 0; } diff --git a/tst/testrandom.g b/tst/testrandom.g index 2407bcf057..5da62128bf 100644 --- a/tst/testrandom.g +++ b/tst/testrandom.g @@ -116,7 +116,7 @@ end;; randomTest := function(collection, randfunc, checkin...) - local sizeone, randchecker; + local sizeone, randchecker, existingPCG32; if Length(checkin) = 0 then checkin := \in; else @@ -145,4 +145,8 @@ randomTest := function(collection, randfunc, checkin...) randchecker(IsGAPRandomSource, GlobalRandomSource, x -> randfunc(GlobalRandomSource, x), randfunc, collection, checkin); + + existingPCG32 := RandomSource(IsPCG32RandomSource); + randchecker(IsPCG32RandomSource, + existingPCG32, x -> randfunc(existingPCG32, x), randfunc, collection, checkin); end;