Pseudo-random numbers/PCG32

Some definitions to help in the explanation

Floor operation

https://en.wikipedia.org/wiki/Floor_and_ceiling_functions

Greatest integer less than or equal to a real number.

Bitwise Logical shift operators (c-inspired)

https://en.wikipedia.org/wiki/Bitwise_operation#Bit_shifts

Binary bits of value shifted left or right, with zero bits shifted in where appropriate.

Examples are shown for 8 bit binary numbers; most significant bit to the left.

<< Logical shift left by given number of bits.

E.g Binary 00110101 << 2 == Binary 11010100

>> Logical shift right by given number of bits.

E.g Binary 00110101 >> 2 == Binary 00001101

^ Bitwise exclusive-or operator

https://en.wikipedia.org/wiki/Exclusive_or

Bitwise comparison for if bits differ

E.g Binary 00110101 ^ Binary 00110011 == Binary 00000110

| Bitwise or operator

https://en.wikipedia.org/wiki/Bitwise_operation#OR

Bitwise comparison gives 1 if any of corresponding bits are 1

E.g Binary 00110101 | Binary 00110011 == Binary 00110111

PCG32 Generator (pseudo-code)

   /* Let u64 denote an unsigned 64 bit integer type. */
   /* Let u32 denote an unsigned 32 bit integer type. */
   
   u64 global constant CONST = 6364136223846793005
   
   class Pcg32
       u64 instance variable state = HEX 853c49e6748fea9b 
       u64 instance variable inc   = HEX da3e39cb94b95bdb    /* Always odd */
   
       method seed(u64 seed_state, u64 seed_sequence)
           state =  0
           inc = (seed_sequence << 1) | 1  # shift and | ensures its odd
           next_int()
           state = state + seed_state
           next_int()
       end method
       
       method next_int()
           u64 old = state
           state = (old * CONST) + inc
           u32 xorshifted = ((old >> 18) ^ old) >> 27
           u32 rot = old >> 59
           u32 answer = (xorshifted >> rot) | (xorshifted << ((-rot) & 31))
           
           return answer
       end method
       
       method next_float():
           return float next_int() / (1 << 32)
       end method
       
   end class

PCG32 Use

   random_gen = instance Pcg_32
   random_gen.seed(42, 54)
   print(random_gen.next_int())   /* 2707161783 */
   print(random_gen.next_int())   /* 2068313097 */
   print(random_gen.next_int())   /* 3122475824 */
   print(random_gen.next_int())   /* 2211639955 */
   print(random_gen.next_int())   /* 3215226955 */

Task

Generate a class/set of functions that generates pseudo-random

numbers as shown above.

Show that the first five integers genrated with the seed 42, 54

are as shown above

Show that for an initial seed of 987654321, 1 the counts of 100_000 repetitions of

   floor(random_gen.next_float() * 5)

Is as follows:

   0: 20049, 1: 20022, 2: 20115, 3: 19809, 4: 20005

Show your output here, on this page.

Delphi

Library: System.SysUtils

Library: Velthuis.BigIntegers

Velthuis.BigIntegers[1] by Rudy Velthuis.

Library: System.Generics.Collections

Translation of: Python

<lang Delphi> program PCG32_test;

{$APPTYPE CONSOLE} uses

 System.SysUtils,
 Velthuis.BigIntegers,
 System.Generics.Collections;

type

 TPCG32 = class
 public
   FState: BigInteger;
   FInc: BigInteger;
   mask64: BigInteger;
   mask32: BigInteger;
   k: BigInteger;
   constructor Create(seedState, seedSequence: BigInteger); overload;
   constructor Create(); overload;
   destructor Destroy; override;
   procedure Seed(seed_state, seed_sequence: BigInteger);
   function NextInt(): BigInteger;
   function NextIntRange(size: Integer): TArray<BigInteger>;
   function NextFloat(): Extended;
 end;

{ TPCG32 }

constructor TPCG32.Create(seedState, seedSequence: BigInteger); begin

 Create();
 Seed(seedState, seedSequence);

end;

constructor TPCG32.Create; begin

 k := '6364136223846793005';
 mask64 := (BigInteger(1) shl 64) - 1;
 mask32 := (BigInteger(1) shl 32) - 1;
 FState := 0;
 FInc := 0;

end;

destructor TPCG32.Destroy; begin

 inherited;

end;

function TPCG32.NextFloat: Extended; begin

 Result := (NextInt.AsExtended / (BigInteger(1) shl 32).AsExtended);

end;

function TPCG32.NextInt(): BigInteger; var

 old, xorshifted, rot, answer: BigInteger;

begin

 old := FState;
 FState := ((old * k) + FInc) and mask64;
 xorshifted := (((old shr 18) xor old) shr 27) and mask32;
 rot := (old shr 59) and mask32;
 answer := (xorshifted shr rot.AsInteger) or (xorshifted shl ((-rot) and
   BigInteger(31)).AsInteger);
 Result := answer and mask32;

end;

function TPCG32.NextIntRange(size: Integer): TArray<BigInteger>; var

 i: Integer;

begin

 SetLength(Result, size);
 if size = 0 then
   exit;

 for i := 0 to size - 1 do
   Result[i] := NextInt;

end;

procedure TPCG32.Seed(seed_state, seed_sequence: BigInteger); begin

 FState := 0;
 FInc := ((seed_sequence shl 1) or 1) and mask64;
 nextint();
 Fstate := (Fstate + seed_state);
 nextint();

end;

var

 PCG32: TPCG32;
 i, key: Integer;
 count: TDictionary<Integer, Integer>;

begin

 PCG32 := TPCG32.Create(42, 54);

 for i := 0 to 4 do
   Writeln(PCG32.NextInt().ToString);

 PCG32.seed(987654321, 1);

 count := TDictionary<Integer, Integer>.Create();

 for i := 1 to 100000 do
 begin
   key := Trunc(PCG32.NextFloat * 5);
   if count.ContainsKey(key) then
     count[key] := count[key] + 1
   else
     count.Add(key, 1);
 end;

 Writeln(#10'The counts for 100,000 repetitions are:');

 for key in count.Keys do
   Writeln(key, ' : ', count[key]);

 count.free;
 PCG32.free;
 Readln;

end. </lang>

Output:

2707161783
2068313097
3122475824
2211639955
3215226955

The counts for 100,000 repetitions are:
3 : 19809
0 : 20049
4 : 20005
2 : 20115
1 : 20022

<lang fsharp> // PCG32. Nigel Galloway: August 13th., 2020 let N=6364136223846793005UL let seed n g=let g=g<<<1|||1UL in (g,(g+n)*N+g) let pcg32=Seq.unfold(fun(n,g)->let rot,xs=uint32(g>>>59),uint32(((g>>>18)^^^g)>>>27) in Some(uint32((xs>>>(int rot))|||(xs<<<(-(int rot)&&&31))),(n,g*N+n))) let pcgFloat n=pcg32 n|>Seq.map(fun n-> (float n)/4294967296.0) </lang>

The Tasks

<lang fsharp> pcg32(seed 42UL 54UL)|>Seq.take 5|>Seq.iter(printfn "%d") </lang>

Output:

<lang fsharp> pcgFloat(seed 987654321UL 1UL)|>Seq.take 100000|>Seq.countBy(fun n->int(n*5.0))|>Seq.iter(printf "%A");printfn "" </lang>

(2, 20115)(3, 19809)(0, 20049)(4, 20005)(1, 20022)

Factor

Translation of: Python

Translation of: Julia

<lang factor>USING: accessors kernel locals math math.statistics prettyprint sequences ;

CONSTANT: mask64 0xffffffffffffffff CONSTANT: mask32 0xffffffff CONSTANT: const 6364136223846793005

TUPLE: pcg32 state inc ;

<pcg32> ( -- pcg32 )

   0x853c49e6748fea9b 0xda3e39cb94b95bdb pcg32 boa ;

next-int ( pcg -- n )

   pcg state>> :> old
   old const * pcg inc>> + mask64 bitand pcg state<<
   old -18 shift old bitxor -27 shift mask32 bitand :> shifted
   old -59 shift mask32 bitand :> r
   shifted r neg shift
   shifted r neg 31 bitand shift bitor mask32 bitand ;

next-float ( pcg -- x ) next-int 1 32 shift /f ;

seed ( pcg st seq -- )

   0x0 pcg state<<
   seq 0x1 shift 1 bitor mask64 bitand pcg inc<<
   pcg next-int drop
   pcg state>> st + pcg state<<
   pcg next-int drop ;

! Task <pcg32> 42 54 [ seed ] keepdd 5 [ dup next-int . ] times

987654321 1 [ seed ] keepdd 100,000 [ dup next-float 5 * >integer ] replicate nip histogram .</lang>

Output:

2707161783
2068313097
3122475824
2211639955
3215226955
H{ { 0 20049 } { 1 20022 } { 2 20115 } { 3 19809 } { 4 20005 } }

Go

Translation of: Python

<lang go>package main

import (

   "fmt"
   "math"

)

const CONST = 6364136223846793005 const mask32 = (1 << 32) - 1

type Pcg32 struct{ state, inc uint64 }

func Pcg32New() *Pcg32 { return &Pcg32{0x853c49e6748fea9b, 0xda3e39cb94b95bdb} }

func (pcg *Pcg32) seed(seedState, seedSequence uint64) {

   pcg.state = 0
   pcg.inc = (seedSequence << 1) | 1
   pcg.nextInt()
   pcg.state = pcg.state + seedState
   pcg.nextInt()

}

func (pcg *Pcg32) nextInt() uint32 {

   old := pcg.state
   pcg.state = old*CONST + pcg.inc
   pcgshifted := uint32((((old >> 18) ^ old) >> 27) & mask32)
   rot := uint32((old >> 59) & mask32)
   return (pcgshifted >> rot) | (pcgshifted << ((-rot) & 31))

}

func (pcg *Pcg32) nextFloat() float64 {

   return float64(pcg.nextInt()) / (1 << 32)

}

func main() {

   randomGen := Pcg32New()
   randomGen.seed(42, 54)
   for i := 0; i < 5; i++ {
       fmt.Println(randomGen.nextInt())
   }

   var counts [5]int
   randomGen.seed(987654321, 1)
   for i := 0; i < 1e5; i++ {
       j := int(math.Floor(randomGen.nextFloat() * 5))
       counts[j]++
   }
   fmt.Println("\nThe counts for 100,000 repetitions are:")
   for i := 0; i < 5; i++ {
       fmt.Printf("  %d : %d\n", i, counts[i])
   }

}</lang>

Output:

2707161783
2068313097
3122475824
2211639955
3215226955

The counts for 100,000 repetitions are:
  0 : 20049
  1 : 20022
  2 : 20115
  3 : 19809
  4 : 20005

Julia

Translation of: Python

<lang julia>const mask64, mask32, CONST = 0xffffffffffffffff, 0xffffffff, UInt(6364136223846793005)

mutable struct PCG32

   state::UInt64
   inc::UInt64
   PCG32(st=0x853c49e6748fea9b, i=0xda3e39cb94b95bdb) = new(st, i)

end

"""return random 32 bit unsigned int""" function next_int!(x::PCG32)

   old = x.state
   x.state = ((old * CONST) + x.inc) & mask64
   xorshifted = (((old >> 18) ⊻ old) >> 27) & mask32
   rot = (old >> 59) & mask32
   return ((xorshifted >> rot) | (xorshifted << ((-rot) & 31))) & mask32

end

"""return random float between 0 and 1""" next_float!(x::PCG32) = next_int!(x) / (1 << 32)

function seed!(x::PCG32, st, seq)

   x.state = 0x0
   x.inc = ((UInt(seq) << 0x1) | 1) & mask64
   next_int!(x)
   x.state = x.state + UInt(st)
   next_int!(x)

end

function testPCG32()

   random_gen = PCG32()
   seed!(random_gen, 42, 54)
   for _ in 1:5
       println(next_int!(random_gen))
   end
   seed!(random_gen, 987654321, 1)
   hist = fill(0, 5)
   for _ in 1:100_000
       hist[Int(floor(next_float!(random_gen) * 5)) + 1] += 1
   end
   println(hist)
   for n in 1:5
       print(n - 1, ": ", hist[n], "  ")
   end

end

testPCG32()

</lang>

Output:

2707161783
2068313097
3122475824
2211639955
3215226955
[20049, 20022, 20115, 19809, 20005]
0: 20049  1: 20022  2: 20115  3: 19809  4: 20005

Phix

Phix proudly does not support the kind of "maths" whereby 255 plus 1 is 0 (or 127+1 is -128).
Phix atoms are limited to 53/64 bits of precision, however (given the above) this task would need 128 bits.
First, for comparison only, this is the usual recommended native approach for this genre of task (different output) <lang Phix>puts(1,"NB: These are not expected to match the task spec!\n") set_rand(42) for i=1 to 5 do

   printf(1,"%d\n",rand(-1))

end for set_rand(987654321) sequence s = repeat(0,5) for i=1 to 100000 do

   s[floor(rnd()*5)+1] += 1

end for ?s</lang>

Output:

NB: These are not expected to match the task spec!
13007222
848581373
2714853861
808614160
2634828316
{20080,19802,19910,20039,20169}

To meet the spec, similar to the Delphi and Wren entries, we resort to using mpfr/gmp, but it is a fair bit longer, and slower, than the above. <lang Phix>include mpfr.e mpz const = mpz_init("6364136223846793005"),

   state = mpz_init("0x853c49e6748fea9b"),
     inc = mpz_init("0xda3e39cb94b95bdb"),   /* Always odd */
     b64 = mpz_init("0x10000000000000000"),  -- (truncate to 64 bits)
     b32 = mpz_init("0x100000000"),          -- (truncate to 32 bits)
     old = mpz_init(),
   xorsh = mpz_init()

procedure seed(integer seed_state, seed_sequence)

   mpz_set_si(inc,seed_sequence*2+1)
   -- as per the talk page:
   -- state := remainder((inc+seed_state)*const+inc,b64)
   mpz_add_ui(state,inc,seed_state)
   mpz_mul(state,state,const)
   mpz_add(state,state,inc)
   mpz_fdiv_r(state, state, b64) -- state := remainder(state,b64)

end procedure

function next_int()

   mpz_set(old,state)          -- old := state
   mpz_set(state,inc)          -- state := inc
   mpz_addmul(state,old,const) -- state += old*const
   mpz_fdiv_r(state, state, b64) -- state := remainder(state,b64) 
   mpz_tdiv_q_2exp(xorsh, old, 18) -- xorsh := trunc(old/2^18)
   mpz_xor(xorsh, xorsh, old)      -- xorsh := xor_bits(xorsh,old)
   mpz_tdiv_q_2exp(xorsh, xorsh, 27) -- xorsh := trunc(xorsh/2^17)
   mpz_fdiv_r(xorsh, xorsh, b32)     -- xorsh := remainder(xorsh,b32) 
   atom xorshifted = mpz_get_atom(xorsh)
   mpz_tdiv_q_2exp(old, old, 59)           -- old := trunc(old/2^59)
   integer rot = mpz_get_integer(old)
   atom answer = xor_bitsu((xorshifted >> rot),(xorshifted << 32-rot))
   return answer

end function

function next_float()

   return next_int() / (1 << 32)

end function

seed(42, 54) for i=1 to 5 do

   printf(1,"%d\n",next_int())

end for seed(987654321,1) sequence r = repeat(0,5) for i=1 to 100000 do

   r[floor(next_float()*5)+1] += 1

end for ?r</lang>

Output:

2707161783
2068313097
3122475824
2211639955
3215226955
{20049,20022,20115,19809,20005}

Python

<lang python>mask64 = (1 << 64) - 1 mask32 = (1 << 32) - 1 CONST = 6364136223846793005

class PCG32():

   def __init__(self, seed_state=None, seed_sequence=None):
       if all(type(x) == int for x in (seed_state, seed_sequence)):
           self.seed(seed_state, seed_sequence)
       else:
           self.state = self.inc = 0
   
   def seed(self, seed_state, seed_sequence):
       self.state = 0
       self.inc = ((seed_sequence << 1) | 1) & mask64
       self.next_int()
       self.state = (self.state + seed_state)
       self.next_int()
       
   def next_int(self):
       "return random 32 bit unsigned int"
       old = self.state
       self.state = ((old * CONST) + self.inc) & mask64
       xorshifted = (((old >> 18) ^ old) >> 27) & mask32
       rot = (old >> 59) & mask32
       answer = (xorshifted >> rot) | (xorshifted << ((-rot) & 31))
       answer = answer &mask32
       
       return answer
   
   def  next_float(self):
       "return random float between 0 and 1"
       return self.next_int() / (1 << 32)

if __name__ == '__main__':

   random_gen = PCG32()
   random_gen.seed(42, 54)
   for i in range(5):
       print(random_gen.next_int())
       
   random_gen.seed(987654321, 1)
   hist = {i:0 for i in range(5)}
   for i in range(100_000):
       hist[int(random_gen.next_float() *5)] += 1
   print(hist)</lang>

Output:

2707161783
2068313097
3122475824
2211639955
3215226955
{0: 20049, 1: 20022, 2: 20115, 3: 19809, 4: 20005}

Raku

Works with: Rakudo version 2020.07

Translation of: Python

Or... at least, it started out that way.

Raku does not have unsigned Integers at this time (Integers are arbitrary sized) so use explicit bit masks during bitwise operations.

<lang perl6>class PCG32 {

   has $!state;
   has $!incr;
   constant mask32 = 2³² - 1;
   constant mask64 = 2⁶⁴ - 1;
   constant const = 6364136223846793005;

   submethod BUILD (
       Int :$seed = 0x853c49e6748fea9b, # default seed
       Int :$incr = 0xda3e39cb94b95bdb  # default increment
     ) {
       $!incr  = $incr +< 1 +| 1 +& mask64;
       $!state = (($!incr + $seed) * const + $!incr) +& mask64;
   }

   method next-int {
       my $shift  = ($!state +> 18 +^ $!state) +> 27 +& mask32;
       my $rotate =  $!state +> 59 +& 31;
       $!state    = ($!state * const + $!incr) +& mask64;
       ($shift +> $rotate) +| ($shift +< (32 - $rotate) +& mask32)
   }

   method next-rat { self.next-int / 2³² }

}

Test next-int with custom seed and increment

say 'Seed: 42, Increment: 54; first five Int values:'; my $rng = PCG32.new( :seed(42), :incr(54) ); .say for $rng.next-int xx 5;

Test next-rat (since these are rational numbers by default)

say "\nSeed: 987654321, Increment: 1; first 1e5 Rat values histogram:"; $rng = PCG32.new( :seed(987654321), :incr(1) ); say ( ($rng.next-rat * 5).floor xx 100_000 ).Bag;

Test next-int with default seed and increment

say "\nSeed: default, Increment: default; first five Int values:"; $rng = PCG32.new; .say for $rng.next-int xx 5;</lang>

Output:

Seed: 42, Increment: 54; first five Int values:
2707161783
2068313097
3122475824
2211639955
3215226955

Seed: 987654321, Increment: 1; first 1e5 Rat values histogram:
Bag(0(20049), 1(20022), 2(20115), 3(19809), 4(20005))

Seed: default, Increment: default; first five Int values:
465482994
3895364073
1746730475
3759121132
2984354868

Wren

Translation of: Python

Library: Wren-big

As Wren doesn't have a 64-bit integer type, we use BigInt instead. <lang ecmascript>import "/big" for BigInt

var Const = BigInt.new("6364136223846793005") var Mask64 = (BigInt.one << 64) - BigInt.one var Mask32 = (BigInt.one << 32) - BigInt.one

class Pcg32 {

   construct new() {
       _state  = BigInt.fromBaseString("853c49e6748fea9b", 16)
       _inc    = BigInt.fromBaseString("da3e39cb94b95bdb", 16)
   }

   seed(seedState, seedSequence) {
       _state = BigInt.zero
       _inc = ((seedSequence << BigInt.one) | BigInt.one) & Mask64
       nextInt
       _state = _state + seedState
       nextInt
   }

   nextInt {
       var old = _state
       _state = (old*Const + _inc) & Mask64
       var xorshifted = (((old >> 18) ^ old) >> 27) & Mask32
       var rot = (old >> 59) & Mask32
       return ((xorshifted >> rot) | (xorshifted << ((-rot) & 31))) & Mask32
   }

   nextFloat { nextInt.toNum / 2.pow(32) }

}

var randomGen = Pcg32.new() randomGen.seed(BigInt.new(42), BigInt.new(54)) for (i in 0..4) System.print(randomGen.nextInt)

var counts = List.filled(5, 0) randomGen.seed(BigInt.new(987654321), BigInt.one) for (i in 1..1e5) {

   var i = (randomGen.nextFloat * 5).floor
   counts[i] = counts[i] + 1

} System.print("\nThe counts for 100,000 repetitions are:") for (i in 0..4) System.print(" %(i) : %(counts[i])")</lang>

Output:

2707161783
2068313097
3122475824
2211639955
3215226955

The counts for 100,000 repetitions are:
  0 : 20049
  1 : 20022
  2 : 20115
  3 : 19809
  4 : 20005

Pseudo-random numbers/PCG32

Delphi

F#

The Functions

The Tasks

Factor

Go

Julia

Phix

Python

Raku

Wren