22
$\begingroup$

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")?

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – ciao Aug 12 '15 at 22:32
  • $\begingroup$ 64-bit, signed, in recent 64-bit versions. 32-bit, signed, for 32-bit versions and all versions before 9. $\endgroup$ – Oleksandr R. Aug 13 '15 at 0:21
  • 2
    $\begingroup$ 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. $\endgroup$ – J. M. will be back soon Aug 13 '15 at 1:55
  • 2
    $\begingroup$ @blochwave I implemented a LibraryLink version of this but I'm starting to lose my confidence in the generation framework, so I got stuck ... $\endgroup$ – Szabolcs Aug 13 '15 at 9:30
  • 1
    $\begingroup$ No reference that I can remember, sorry. But you can verify it using Developer`MachineIntegerQ. $\endgroup$ – Oleksandr R. Aug 13 '15 at 10:19
20
$\begingroup$

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

histograms for the PCG generator at various bit widths

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]]]

"noise sphere" for PCG generator

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – J. M. will be back soon Aug 17 '15 at 14:55
  • 1
    $\begingroup$ This is really cool, and really thorough - thank you! $\endgroup$ – dr.blochwave Aug 17 '15 at 15:33
16
$\begingroup$

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

enter image description here

$\endgroup$
  • 1
    $\begingroup$ 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? $\endgroup$ – J. M. will be back soon Aug 14 '15 at 3:52
  • $\begingroup$ @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. $\endgroup$ – Sasha Aug 14 '15 at 19:17
  • $\begingroup$ Amazing - this looks good, thank you! $\endgroup$ – dr.blochwave Aug 14 '15 at 20:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.