# Implementing the PCG random number generator

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.

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 RandomInitializeGenerator[gsym, opts] (as per "Defining your own generator")?

• So, the question is "...is it then possible to plug it into the random framework as our own custom RNG..."? Then the answer is yes.
– ciao
Aug 12, 2015 at 22:32
• 64-bit, signed, in recent 64-bit versions. 32-bit, signed, for 32-bit versions and all versions before 9. Aug 13, 2015 at 0:21
• This looks interesting. I've some experience with implementing new PRNG algorithms with the provided framework; I'll give this a whirl on the weekend. Aug 13, 2015 at 1:55
• @blochwave I implemented a LibraryLink version of this but I'm starting to lose my confidence in the generation framework, so I got stuck ... Aug 13, 2015 at 9:30
• No reference that I can remember, sorry. But you can verify it using DeveloperMachineIntegerQ. Aug 13, 2015 at 10:19

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 /: RandomInitializeGenerator[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[]]]; inc = "Increment" /. flops; If[inc === Automatic, If[1 <= idx <= 5, inc =$PCGDefaultIncrements[[idx]],
Message[PermutedCongruential::ni, bw]; bw = 32;
inc = $PCGDefaultIncrements[]], 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. • From the looks of it, the only thing preventing the extension of PCG to 128 bits and higher is the lack of good (M)LCGs that are 256 bits or higher. Unfortunately, L'Ecuyer's listings only go up to 128 bits. Aug 17, 2015 at 14:55 • This is really cool, and really thorough - thank you! Aug 17, 2015 at 15:33 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 /: RandomInitializeGenerator[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:= BlockRandom[
SeedRandom[12324235,
Method -> {PermutedCongruential, "InitSeq" -> 32736465}];
{RandomInteger[{1, 26}, 3], RandomReal[1, 3]}]

Out= {{21, 9, 15}, {0.175217, 0.750743, 0.914359}}


The results appear to be good:

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

Out= 0.925037 • I think it's related to the point I raised in a comment: how are you supposed to "correctly" (whatever that means) rotate bits in Mathematica? Aug 14, 2015 at 3:52
• @J.M. The issue was incorrect "BitWidth". I used 64, and it gives 32, which explained the atom at 0 I was seeing. The other issue was that I needed to truncate bit-expandable Mathematica integers to the proper widths as in C-code. I do this by BitAnd-ing with a mask now. Aug 14, 2015 at 19:17
• Amazing - this looks good, thank you! Aug 14, 2015 at 20:35