Implementing the PCG random number generator

Question

PCG ("permuted congruential generator") is a recently-introduced random number generator.

PCG is a family of simple fast space-efficient statistically good algorithms for random number generation. Unlike many general-purpose RNGs, they are also hard to predict.

Source: PCG, A Family of Better Random Number Generators

There's a minimal 9-line implementation in C, along with C, C++, Haskell and also Python implementations.

How can we implement the PCG directly in Mathematica?

I've translated the 9-line C into Mathematica, but I'm now very confused thanks to Mathematica's internal representation of integers (32-bit? 64-bit? N-bit?) so I'm not sure I've done it correctly, let alone optimally.

pcgRandomR[state_, inc_] := 
 Module[{oldstate, newstate, xorshifted, rot},
    oldstate = state;
    newstate = oldstate*6364136223846793005 + BitOr[inc, 1];
    xorshifted = BitShiftRight[BitXor[BitShiftRight[oldstate, 18], oldstate], 27];
    rot = BitShiftRight[oldstate, 59];
    {BitOr[BitShiftRight[xorshifted, rot], 
        BitShiftLeft[xorshifted, BitAnd[-rot, 31]]], newstate, inc}
  ]

This is my translation of the seeding function found here on PCG.SE. Again, not really sure about correctness due to the representation of integers.

pcgRandomRSeed[initstate_, initseq_] := Module[{state, inc, rand},
     state = 0;
     inc = BitOr[BitShiftLeft[initseq, 1], 1];
     {rand, state, inc} = pcgRandomR[state, inc];
     state += initstate;
     {rand, state, inc} = pcgRandomR[state, inc];
     {state, inc}
  ]

Having sorted out those problems in defining the PCG, is it then possible to plug it into the random framework as our own custom RNG using Random`InitializeGenerator[gsym, opts] (as per "Defining your own generator")?

Answer 1

After some amount of effort, I managed to come up with an implementation of O'Neill's "XSH-RR" family of permuted congruential generators. The following covers the 8-, 16-, 32-, and 64-bit generators, and the mcg, oneseq, and setseq variants. (I'll leave the modification to handle the unique variant as an exercise for the interested reader.) A similar approach can be done with the other PCG families.

First, here is a routine to rotate bits to the right:

BitRotateRight[n_Integer, r_Integer?NonNegative, bits : (_Integer?Positive) : 32] /;
               BitAnd[bits, bits - 1] == 0 := BitShiftRight[n, r] + 
               BitShiftLeft[BitAnd[n, BitShiftLeft[1, r] - 1], bits - r]

(This effectively generalizes Daniel's code here.)

Here is a set of default multipliers and increments for the underlying linear congruential generators. These and the following were adapted from O'Neill's code here:

$PCGDefaultMultipliers = {141, 12829, 747796405, 6364136223846793005, 
                          47026247687942121848144207491837523525};

$PCGDefaultIncrements = {77, 47989, 2891336453, 1442695040888963407, 
                         117397592171526113268558934119004209487};

A few error messages:

PermutedCongruential::bw = "The value of the option \"BitWidth\"\[Rule]`1`
should be a power of two or Automatic.";
PermutedCongruential::nm = "No built-in multiplier available for
\"BitWidth\"\[Rule]`1`. Using \"BitWidth\"\[Rule]32 instead.";
PermutedCongruential::ni = "No built-in increment available for
\"BitWidth\"\[Rule]`1`. Using \"BitWidth\"\[Rule]32 instead.";

Some options:

Options[PermutedCongruential] =
  {"BitWidth" -> Automatic, "Increment" -> Automatic, "Multiplier" -> Automatic};

Initialize the generator:

PermutedCongruential /: Random`InitializeGenerator[PermutedCongruential,
                                                   opts___] := 
Module[{flops = Flatten[{opts, Options[PermutedCongruential]}],
        bw, idx, inc, mul},

       bw = "BitWidth" /. flops;
       If[bw === Automatic, bw = 32, 
          If[! (IntegerQ[bw] && Positive[bw] && BitAnd[bw, bw - 1] == 0), 
             Message[PermutedCongruential::bw, bw]; Throw[$Failed]]];

       idx = BitLength[bw] - 2;

       mul = "Multiplier" /. flops;
       If[mul === Automatic, 
          If[1 <= idx <= 5, mul = $PCGDefaultMultipliers[[idx]], 
             Message[PermutedCongruential::nm, bw]; bw = 32;
             mul = $PCGDefaultMultipliers[[4]]]];

       inc = "Increment" /. flops;
       If[inc === Automatic, 
          If[1 <= idx <= 5, inc = $PCGDefaultIncrements[[idx]], 
             Message[PermutedCongruential::ni, bw]; bw = 32;
             inc = $PCGDefaultIncrements[[4]]], 
          inc = BitOr[BitShiftLeft[inc], 1]];

       idx = Min[Max[1, idx], 5];
       PermutedCongruential[bw, mul, inc, {BitShiftRight[bw + idx + 1],
                            bw - idx - 1, 2 bw - idx - 1}, 0]]

We define PCG as a bit generator:

PermutedCongruential[___]["GeneratesBitsQ"] := True
PermutedCongruential[bw_, rest__]["BitWidth"] := bw

This handles seeding:

PermutedCongruential[bw_, mul_, inc_, const_, state_]["SeedGenerator"[seed_]] :=
    Module[{ss = seed + BitAnd[inc, BitShiftLeft[1, 2 bw] - 1]},
           PermutedCongruential[bw, mul, inc, const, 
                                BitAnd[mul ss + inc, BitShiftLeft[1, 2 bw] - 1]]]

Finally, here is the bit generator:

PermutedCongruential[bw_, mul_, inc_, const_, state_]["GenerateBits"[bits_]] :=
    Module[{bot, ns, res, rot, xsh},
           {xsh, bot, rot} = const;
           ns = BitAnd[mul state + inc, BitShiftLeft[1, 2 bw] - 1];
           res = BitShiftRight[BitXor[state, BitShiftRight[state, xsh]], bot];
           res = BitRotateRight[BitAnd[res, BitShiftLeft[1, bw] - 1],
                                BitShiftRight[state, rot], bw];
           {res, PermutedCongruential[bw, mul, inc, const, ns]}]

---

As a test of the generator, let's reproduce some sample values given here:

BlockRandom[SeedRandom[42, Method -> {PermutedCongruential, "Increment" -> 54}];
            IntegerString[RandomInteger[{0, BitShiftLeft[1, 32] - 1}, 6], 16]]
   {"a15c02b7", "7b47f409", "ba1d3330", "83d2f293", "bfa4784b", "cbed606e"}

Test uniformity at various bit widths:

Table[DistributionFitTest[BlockRandom[SeedRandom[42,
      Method -> {PermutedCongruential, "BitWidth" -> BitShiftLeft[1, k]}];
      RandomReal[1, 1*^5]], UniformDistribution[]], {k, 3, 6}]
   {0.988711, 0.951478, 0.228108, 0.181171}

Table[Histogram[BlockRandom[SeedRandom[42, 
      Method -> {PermutedCongruential, "BitWidth" -> BitShiftLeft[1, k]}];
      RandomReal[1, 1*^5]], Automatic, "PDF", 
      PlotLabel -> StringForm["``-bit PCG", BitShiftLeft[1, k]]],
      {k, 3, 6}] // Partition[#, 2] & // GraphicsGrid

At this point, I'd like to take this opportunity to present another method to visualize the randomness of a sequence. This is Pickover's "noise sphere" visualization:

noiseSphere[vals_?VectorQ, opts___] :=
     Graphics3D[{RGBColor[##], 
                 Sphere[Sqrt[#3] Append[Through[{Cos, Sin}[2 π #1]] Sin[π #2], 
                        Cos[π #2]], 1/50]} & @@@ Partition[vals, 3, 1],
                opts, Boxed -> False, Lighting -> "Neutral"]

BlockRandom[SeedRandom[42, Method -> {PermutedCongruential, "Increment" -> 54}];
            noiseSphere[RandomReal[1, 1*^4]]]

The lack of structure and uniform filling of the sphere visually demonstrate the (apparent) lack of correlation in the sequence generated by PCG.

Answer 2

I just followed through with the tutorial "Defining Your Own Generator".

Start with provided, a little tweaked functions. The key trick is to ensure that the bit-expandable Mathematica integers are of the size of relative machine unsigned integers. I use BitAnd with mask to accomplish that:

pcgRandomR[state_, inc_] := 
 Module[{ newstate, xorshifted, rot, mask32=2^32-1,mask64=2^64-1},
    newstate = BitAnd[state*6364136223846793005 + BitOr[inc, 1], mask64];
    xorshifted = BitAnd[BitShiftRight[
        BitXor[BitShiftRight[state, 18], state], 27], mask32];
    rot = Mod[BitShiftRight[state, 59], 32];
    {BitAnd[BitOr[
        BitShiftRight[xorshifted, rot], 
        BitShiftLeft[xorshifted, BitAnd[32 - rot, 31]
     ]], mask32], newstate, inc}
  ]

pcgRandomRSeed[initstate_, initseq_] := Module[{state, inc, rand},
     state = 0;
     inc = Mod[BitOr[BitShiftLeft[initseq, 1], 1], 2^64];
     {rand, state, inc} = pcgRandomR[state, inc];
     state = Mod[state + initstate, 2^64];
     {rand, state, inc} = pcgRandomR[state, inc];
     {state, inc}
  ]

Plug into the framework:

Options[PermutedCongruential] = {"InitState" -> Automatic, 
  "InitSeq" -> 7050702485517258437}

PermutedCongruential /:
 Random`InitializeGenerator[PermutedCongruential, opts___] :=
 Module[{initState, initSeq},
  initState = 
   Replace[OptionValue[PermutedCongruential, {opts}, "InitState"], 
    Automatic :> RandomInteger[{1, 2^64}]];
  initSeq = 
   Replace[OptionValue[PermutedCongruential, {opts}, "InitSeq"], 
    Automatic :> RandomInteger[{1, 2^64}]];
  If[! IntegerQ[initState], Throw[$Failed]];
  If[! IntegerQ[initSeq], Throw[$Failed]];
  PermutedCongruential @@ pcgRandomRSeed[initState, initSeq]
  ]

The generated random integers are 32-bits, hence the bit-width of the generator must be 32, rather 64, as was the case in my previous unsuccessful attempt at it:

PermutedCongruential[state_, inc_]["GeneratesBitsQ"] := True;
PermutedCongruential[state_, inc_]["BitWidth"] = 32;
PermutedCongruential[state_, inc_]["SeedGenerator"[seed_]] :=
 PermutedCongruential[Mod[state seed, 2^64], inc]

PermutedCongruential[state_, inc_]["GenerateBits"[bits_]] := 
    {#1, PermutedCongruential[##2]} & @@ pcgRandomR[state, inc]

It can now be used:

In[78]:= BlockRandom[
 SeedRandom[12324235, 
  Method -> {PermutedCongruential, "InitSeq" -> 32736465}];
 {RandomInteger[{1, 26}, 3], RandomReal[1, 3]}]

Out[78]= {{21, 9, 15}, {0.175217, 0.750743, 0.914359}}

The results appear to be good:

In[79]:= DistributionFitTest[BlockRandom[
  SeedRandom[12324235, 
   Method -> {PermutedCongruential, "InitSeq" -> 32736465}]; 
  RandomReal[1, 10^5]], UniformDistribution[]]

Out[79]= 0.925037