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;