Add PCG32 random number generator

lib/random.gd

@@ -146,11 +146,12 @@ DeclareOperation( "Init", [IsRandomSource, IsObject] );
 ##  <Filt Name="IsMersenneTwister" Arg='rs' Type='Category'/>
 ##  <Filt Name="IsGAPRandomSource" Arg='rs' Type='Category'/>
 ##  <Filt Name="IsGlobalRandomSource" Arg='rs' Type='Category'/>
+##  <Filt Name="IsPCG32RandomSource" Arg='rs' Type='Category'/>
 ##  <Var Name="GlobalMersenneTwister"/>
 ##  <Var Name="GlobalRandomSource"/>
 ##
 ##  <Description>
-##  Currently, the &GAP; library provides three types of random sources,
+##  Currently, the &GAP; library provides four types of random sources,
 ##  distinguished by the three listed categories.
 ##  <P/>
 ##  <Ref  Var="IsMersenneTwister"/>  are random  sources  which  use a  fast
@@ -160,8 +161,14 @@ DeclareOperation( "Init", [IsRandomSource, IsObject] );
 ##  and the origin of the code used in the &GAP; kernel, see:
 ##  <URL>http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html</URL>.
 ##  <P/>
-##  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.
+##  <P/>
+##  <Ref Var="IsPCG32RandomSource"/> 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 <URL>http://www.pcg-random.org/</URL>.
 ##  <P/>
 ##  There is also a predefined global random source 
 ##  <Ref Var="GlobalMersenneTwister"/> 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" );

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.

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;
 }
 

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;