Pseudo-random numbers/PCG32: Difference between revisions
Thundergnat (talk | contribs) (→{{header|Raku}}: More commentary, test default seed and increment too) |
(realize in F#) |
||
Line 234: | Line 234: | ||
</pre> |
</pre> |
||
=={{header|F_Sharp|F#}}== |
|||
===The Functions=== |
|||
<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> |
|||
{{out}} |
|||
<pre> |
|||
2707161783 |
|||
2068313097 |
|||
3122475824 |
|||
2211639955 |
|||
3215226955 |
|||
</pre> |
|||
<lang fsharp> |
|||
pcgFloat(seed 987654321UL 1UL)|>Seq.take 100000|>Seq.countBy(fun n->int(n*5.0))|>Seq.iter(printf "%A");printfn "" |
|||
</lang> |
|||
<pre> |
|||
(2, 20115)(3, 19809)(0, 20049)(4, 20005)(1, 20022) |
|||
</pre> |
|||
=={{header|Factor}}== |
=={{header|Factor}}== |
||
{{trans|Python}} |
{{trans|Python}} |
Revision as of 14:15, 13 August 2020
- 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.
- https://en.wikipedia.org/wiki/Floor_and_ceiling_functions
- 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.
- https://en.wikipedia.org/wiki/Bitwise_operation#Bit_shifts
- << Logical shift left by given number of bits.
- E.g Binary 00110101 << 2 == Binary 11010100
- << Logical shift left by given number of bits.
- >> Logical shift right by given number of bits.
- E.g Binary 00110101 >> 2 == Binary 00001101
- >> Logical shift right by given number of bits.
- ^ 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
Velthuis.BigIntegers[1] by Rudy Velthuis.
<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
F#
The Functions
<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:
2707161783 2068313097 3122475824 2211639955 3215226955
<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
<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
<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
<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
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
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
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