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 would need 128 bits.
First, for comparison only, this is the usual recommended native approach (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()

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(xorsh/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

procedure seed(integer seed_state, seed_sequence) -- mpz_set_si(state,0)

   mpz_set_si(inc,seed_sequence*2+1)

-- {} = next_int() -- mpz_add_ui(state,state,seed_state) -- {} = next_int()

   -- instead, 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_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