Pseudo-random numbers/Xorshift star

From Rosetta Code
Revision as of 21:15, 14 August 2020 by PureFox (talk | contribs) (→‎{{header|Go}}: Rely on implicit conversion from uint64 to uint32 rather than use a mask.)
Pseudo-random numbers/Xorshift star is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
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:
3540625527
2750739987
4037983143
1993361440
3809424708

<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

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)

}

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;

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

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

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

}

  1. Test next-int

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


  1. 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;


  1. 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
Retrieved from "https://rosettacode.org/wiki/Pseudo-random_numbers/Xorshift_star?oldid=301827"