Pseudo-random numbers/Xorshift star: Difference between revisions

Content added Content deleted

Inline

Revision as of 16:42, 14 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.

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

Xorshift_star Generator (pseudo-code)

   /* Let u64 denote an unsigned 64 bit integer type. */
   /* Let u32 denote an unsigned 32 bit integer type. */

   class Xorshift_star
       u64 state       /* Must be seeded to non-zero initial value */
       u64 const = HEX '2545F4914F6CDD1D'

       method seed(u64 num):
           state =  num
       end method
       
       method next_int():
           u64 x = state
           x = x ^ (x >> 12)
           x = x ^ (x << 25)
           x = x ^ (x >> 27)
           state = x
           u32 answer = ((x * const) >> 32)
           
           return answer
       end method
       
       method next_float():
           return float next_int() / (1 << 32)
       end method
       
   end class

Xorshift use

   random_gen = instance Xorshift_star
   random_gen.seed(1234567)
   print(random_gen.next_int())   /* 3540625527 */
   print(random_gen.next_int())   /* 2750739987 */
   print(random_gen.next_int())   /* 4037983143 */
   print(random_gen.next_int())   /* 1993361440 */
   print(random_gen.next_int())   /* 3809424708 */

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 1234567

are as shown above

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

   floor(random_gen.next_float() * 5)

Is as follows:

   0: 20103, 1: 19922, 2: 19937, 3: 20031, 4: 20007

Show your output here, on this page.

F#

The Functions

<lang fsharp> // Xorshift star. Nigel Galloway: August 14th., 2020 let fN=(fun(n:uint64)->n^^^(n>>>12))>>(fun n->n^^^(n<<<25))>>(fun n->n^^^(n>>>27)) let Xstar32=Seq.unfold(fun n->let n=fN n in Some(uint32((n*0x2545F4914F6CDD1DUL)>>>32),n)) let XstarF n=Xstar32 n|>Seq.map(fun n->(float n)/4294967296.0) </lang>

The Tasks

<lang fsharp> Xstar32 1234567UL|>Seq.take 5|>Seq.iter(printfn "%d") </lang>

Output:

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

Output:

(4, 20007)(2, 19937)(3, 20031)(0, 20103)(1, 19922)

Factor

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

CONSTANT: mask64 $[ 1 64 shift 1 - ] CONSTANT: mask32 $[ 1 32 shift 1 - ] CONSTANT: const 0x2545F4914F6CDD1D

! Restrict seed value to positive integers. PREDICATE: positive < integer 0 > ; ERROR: seed-nonpositive seed ;

TUPLE: xorshift* { state positive initial: 1 } ;

<xorshift*> ( seed -- xorshift* )

   dup positive? [ seed-nonpositive ] unless
   mask64 bitand xorshift* boa ;

twiddle ( m n -- n ) dupd shift bitxor mask64 bitand ;

next-int ( obj -- n )

   dup state>> -12 twiddle 25 twiddle -27 twiddle tuck swap
   state<< const * mask64 bitand -32 shift mask32 bitand ;

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

! ---=== Task ===--- 1234567 <xorshift*> 5 [ dup next-int . ] times

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

Output:

3540625527
2750739987
4037983143
1993361440
3809424708
H{ { 0 20103 } { 1 19922 } { 2 19937 } { 3 20031 } { 4 20007 } }

Go

Translation of: Python

<lang go>package main

import (

   "fmt"
   "math"

)

const CONST = 0x2545F4914F6CDD1D const mask32 = (1 << 32) - 1

type XorshiftStar struct{ state uint64 }

func XorshiftStarNew(state uint64) *XorshiftStar { return &XorshiftStar{state} }

func (xor *XorshiftStar) seed(state uint64) { xor.state = state }

func (xor *XorshiftStar) nextInt() uint32 {

   x := xor.state
   x = x ^ (x >> 12)
   x = x ^ (x << 25)
   x = x ^ (x >> 27)
   xor.state = x
   return uint32((x * CONST) >> 32 & mask32)

}

func (xor *XorshiftStar) nextFloat() float64 {

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

}

func main() {

   randomGen := XorshiftStarNew(1234567)
   for i := 0; i < 5; i++ {
       fmt.Println(randomGen.nextInt())
   }

   var counts [5]int
   randomGen.seed(987654321)
   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:

3540625527
2750739987
4037983143
1993361440
3809424708

The counts for 100,000 repetitions are:
  0 : 20103
  1 : 19922
  2 : 19937
  3 : 20031
  4 : 20007

Julia

Translation of: Python

<lang julia>const mask64 = (0x1 << 64) - 1 const mask32 = (0x1 << 32) - 1 const CONST = 0x2545F4914F6CDD1D

mutable struct XorShiftStar state::UInt end

XorShiftStar(_seed=0x0) = XorShiftStar(UInt(_seed) & mask64)

seed(x::XorShiftStar, num) = begin x.state = num & mask64 end

"""return random int between 0 and 2**32""" function next_int(x::XorShiftStar)

   x.state = (x.state ⊻ (x.state >> 12)) & mask64
   x.state = (x.state ⊻ (x.state << 25)) & mask64
   x.state = (x.state ⊻ (x.state >> 27)) & mask64
   return (((x.state * CONST) & mask64) >> 32) & mask32

end

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

function testXorShiftStar()

   random_gen = XorShiftStar()
   seed(random_gen, 1234567)
   for i in 1:5
       println(next_int(random_gen))
   end
   seed(random_gen, 987654321)
   hist = fill(0, 5)
   for _ in 1:100_000
       hist[Int(floor(next_float(random_gen) * 5)) + 1] += 1
   end
   foreach(n -> print(n - 1, ": ", hist[n], "  "), 1:5)

end

testXorShiftStar()

</lang>

Output:

3540625527
2750739987
4037983143
1993361440
3809424708
0: 20103  1: 19922  2: 19937  3: 20031  4: 20007

Phix

As per Pseudo-random_numbers/PCG32#Phix, resorting to mpfr/gmp <lang Phix>include mpfr.e mpz const = mpz_init("0x2545F4914F6CDD1D"),

   state = mpz_init(),
     b64 = mpz_init("0x10000000000000000"),  -- (truncate to 64 bits)
     b32 = mpz_init("0x100000000"),          -- (truncate to 32 bits)
     tmp = mpz_init(),
       x = mpz_init()

procedure seed(integer num)

   mpz_set_si(state,num)

end procedure

function next_int()

   mpz_set(x,state)            -- x := state
   mpz_tdiv_q_2exp(tmp, x, 12) -- tmp := trunc(x/2^12)
   mpz_xor(x, x, tmp)          -- x := xor_bits(x,tmp)
   mpz_mul_2exp(tmp, x, 25)    -- tmp := x * 2^25.
   mpz_xor(x, x, tmp)          -- x := xor_bits(x,tmp)
   mpz_fdiv_r(x, x, b64)       -- x := remainder(x,b64) 
   mpz_tdiv_q_2exp(tmp, x, 27) -- tmp := trunc(x/2^27)
   mpz_xor(x, x, tmp)          -- x := xor_bits(x,tmp)
   mpz_fdiv_r(state, x, b64)   -- state := remainder(x,b64) 
   mpz_mul(x,x,const)          -- x *= const
   mpz_tdiv_q_2exp(x, x, 32)   -- x := trunc(x/2^32)
   mpz_fdiv_r(x, x, b32)       -- x := remainder(x,b32) 
   return mpz_get_atom(x)

end function

function next_float()

   return next_int() / (1 << 32)

end function

seed(1234567) for i=1 to 5 do

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

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

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

end for ?r</lang>

Output:

3540625527
2750739987
4037983143
1993361440
3809424708
{20103,19922,19937,20031,20007}

Python

<lang python>mask64 = (1 << 64) - 1 mask32 = (1 << 32) - 1 const = 0x2545F4914F6CDD1D

class Xorshift_star():

   def __init__(self, seed=0):
       self.state = seed & mask64

   def seed(self, num):
       self.state =  num & mask64
   
   def next_int(self):
       "return random int between 0 and 2**32"
       x = self.state
       x = (x ^ (x >> 12)) & mask64
       x = (x ^ (x << 25)) & mask64
       x = (x ^ (x >> 27)) & mask64
       self.state = x
       answer = (((x * const) & mask64) >> 32) & 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 = Xorshift_star()
   random_gen.seed(1234567)
   for i in range(5):
       print(random_gen.next_int())
       
   random_gen.seed(987654321)
   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:

3540625527
2750739987
4037983143
1993361440
3809424708
{0: 20103, 1: 19922, 2: 19937, 3: 20031, 4: 20007}

Raku

Works with: Rakudo version 2020.07

Translation of: Python

Raku does not have unsigned Integers at this time (Integers are arbitrary sized) so use explicit bit masks during bitwise operations. All constants are encapsulated inside the class.

<lang perl6>class Xorshift-star {

   has $!state;

   constant mask64 = 2⁶⁴ - 1;
   constant const = 0x2545F4914F6CDD1D;

   submethod BUILD ( Int :$seed where * > 0 = 1 ) { $!state = $seed +& mask64 }

   method next-int {
       $!state +^= $!state +> 12;
       $!state +^= $!state +< 25 +& mask64;
       $!state +^= $!state +> 27;
       ($!state * const) +> 32 +& (2³² - 1)
   }

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

}

Test next-int

say 'Seed: 1234567; first five Int values'; my $rng = Xorshift-star.new( :seed(1234567) ); .say for $rng.next-int xx 5;

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

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

Test with default seed

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

Output:

Seed: 1234567; first five Int values
3540625527
2750739987
4037983143
1993361440
3809424708

Seed: 987654321; first 1e5 Rat values histogram
Bag(0(20103), 1(19922), 2(19937), 3(20031), 4(20007))

Seed: default; first five Int values
1206177355
2882512552
3117485455
1303648416
241277360

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.fromBaseString("2545F4914F6CDD1D", 16) var Mask64 = (BigInt.one << 64) - BigInt.one var Mask32 = (BigInt.one << 32) - BigInt.one

class XorshiftStar {

   construct new(state) {
       _state  = state & Mask64
   }

   seed(num) { _state = num & Mask64}

   nextInt {
       var x = _state 
       x = (x ^ (x >> 12)) & Mask64
       x = (x ^ (x << 25)) & Mask64
       x = (x ^ (x >> 27)) & Mask64
       _state = x
       return (((x * Const) & Mask64) >> 32) & Mask32
   }

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

}

var randomGen = XorshiftStar.new(BigInt.new(1234567)) for (i in 0..4) System.print(randomGen.nextInt)

var counts = List.filled(5, 0) randomGen.seed(BigInt.new(987654321)) 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:

3540625527
2750739987
4037983143
1993361440
3809424708

The counts for 100,000 repetitions are:
  0 : 20103
  1 : 19922
  2 : 19937
  3 : 20031
  4 : 20007