I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

Pseudo-random numbers/Xorshift star

From Rosetta Code
Task
Pseudo-random numbers/Xorshift star
You are encouraged to solve this task according to the task description, using any language you may know.
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.


ALGOL 68[edit]

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32

Will generate a runtime error if state is not initialised before use.

BEGIN # generate some pseudo random numbers using Xorshift star #
# note that although LONG INT is 64 bits in Algol 68G, LONG BITS is longer than 64 bits #
LONG BITS state;
LONG INT const = ABS LONG 16r2545f4914f6cdd1d;
LONG INT one shl 32 = ABS ( LONG 16r1 SHL 32 );
# sets the state to the specified seed value #
PROC seed = ( LONG INT num )VOID: state := BIN num;
# XOR and assign convenience operator #
PRIO XORAB = 1;
OP XORAB = ( REF LONG BITS x, LONG BITS v )REF LONG BITS: x := ( x XOR v ) AND LONG 16rffffffffffffffff;
# gets the next pseudo random integer #
PROC next int = LONG INT:
BEGIN
LONG BITS x := state;
x XORAB ( x SHR 12 );
x XORAB ( x SHL 25 );
x XORAB ( x SHR 27 );
state := x;
SHORTEN ABS ( 16rffffffff AND BIN ( ABS x * LENG const ) SHR 32 )
END # next int # ;
# gets the next pseudo random real #
PROC next float = LONG REAL: next int / one shl 32;
BEGIN # task test cases #
seed( 1234567 );
print( ( whole( next int, 0 ), newline ) ); # 3540625527 #
print( ( whole( next int, 0 ), newline ) ); # 2750739987 #
print( ( whole( next int, 0 ), newline ) ); # 4037983143 #
print( ( whole( next int, 0 ), newline ) ); # 1993361440 #
print( ( whole( next int, 0 ), newline ) ); # 3809424708 #
# count the number of occurances of 0..4 in a sequence of pseudo random reals scaled to be in [0..5) #
seed( 987654321 );
[ 0 : 4 ]INT counts; FOR i FROM LWB counts TO UPB counts DO counts[ i ] := 0 OD;
TO 100 000 DO counts[ SHORTEN ENTIER ( next float * 5 ) ] +:= 1 OD;
FOR i FROM LWB counts TO UPB counts DO
print( ( whole( i, -2 ), ": ", whole( counts[ i ], -6 ) ) )
OD;
print( ( newline ) )
END
END
Output:
3540625527
2750739987
4037983143
1993361440
3809424708
 0:  20103 1:  19922 2:  19937 3:  20031 4:  20007

F#[edit]

The Functions[edit]

 
// 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)
 

The Tasks[edit]

 
Xstar32 1234567UL|>Seq.take 5|>Seq.iter(printfn "%d")
 
Output:
3540625527
2750739987
4037983143
1993361440
3809424708
 
XstarF 987654321UL|>Seq.take 100000|>Seq.countBy(fun n->int(n*5.0))|>Seq.iter(printf "%A");printfn ""
 
Output:
(4, 20007)(2, 19937)(3, 20031)(0, 20103)(1, 19922)

Factor[edit]

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 .
Output:
3540625527
2750739987
4037983143
1993361440
3809424708
H{ { 0 20103 } { 1 19922 } { 2 19937 } { 3 20031 } { 4 20007 } }

Go[edit]

Translation of: Python
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])
}
}
Output:
3540625527
2750739987
4037983143
1993361440
3809424708

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

Julia[edit]

Translation of: Python
const mask32 = (0x1 << 32) - 1
const CONST = 0x2545F4914F6CDD1D
 
mutable struct XorShiftStar
state::UInt64
end
 
XorShiftStar(_seed=0x0) = XorShiftStar(UInt(_seed))
 
seed(x::XorShiftStar, num) = begin x.state = UInt64(num) end
 
"""return random int between 0 and 2**32"""
function next_int(x::XorShiftStar)
x.state = x.state ⊻ (x.state >> 12)
x.state = x.state ⊻ (x.state << 25)
x.state = x.state ⊻ (x.state >> 27)
return ((x.state * CONST) >> 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()
 
Output:
3540625527
2750739987
4037983143
1993361440
3809424708
0: 20103  1: 19922  2: 19937  3: 20031  4: 20007

Phix[edit]

As per Pseudo-random_numbers/PCG32#Phix, resorting to mpfr/gmp

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
Output:
3540625527
2750739987
4037983143
1993361440
3809424708
{20103,19922,19937,20031,20007}

Python[edit]

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)
Output:
3540625527
2750739987
4037983143
1993361440
3809424708
{0: 20103, 1: 19922, 2: 19937, 3: 20031, 4: 20007}

Raku[edit]

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.

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³² }
}
 
# 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;
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


REXX[edit]

/*REXX program  generates   pseudo─random numbers   using the  XOR─shift─star  method.  */
numeric digits 200 /*ensure enough decimal digs for mult. */
parse arg n reps pick seed1 seed2 . /*obtain optional arguments from the CL*/
if n=='' | n=="," then n= 5 /*Not specified? Then use the default.*/
if reps=='' | reps=="," then reps= 100000 /* " " " " " " */
if pick=='' | pick=="," then pick= 5 /* " " " " " " */
if seed1=='' | seed1=="," then seed1= 1234567 /* " " " " " " */
if seed2=='' | seed2=="," then seed2= 987654321 /* " " " " " " */
const= x2d(2545f4914f6cdd1d) /*initialize the constant to be used. */
o.12= copies(0, 12) /*construct 12 bits of zeros. */
o.25= copies(0, 25) /* " 25 " " " */
o.27= copies(0, 27) /* " 27 " " " */
w= max(3, length(n) ) /*for aligning the left side of output.*/
state= seed1 /* " the state to seed #1. */
do j=1 for n
if j==1 then do; say center('n', w) " pseudo─random number"
say copies('═', w) " ══════════════════════════"
end
say right(j':', w)" " right(commas(next()), 18) /*display a random number*/
end /*j*/
say
if reps==0 then exit 0 /*stick a fork in it, we're all done. */
say center('#', w) " count of pseudo─random #"
say copies('═', w) " ══════════════════════════"
state= seed2 /* " the state to seed #2. */
@.= 0; div= pick / 2**32 /*convert division to inverse multiply.*/
do k=1 for reps
parse value next()*div with _ '.' /*get random #, floor of a "division". */
@._= @._ + 1 /*bump the counter for this random num.*/
end /*k*/
 
do #=0 for pick
say right(#':', w)" " right(commas(@.#), 14) /*show count of a random num.*/
end /*#*/
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
b2d: parse arg ?; return x2d( b2x(?) ) /*convert bin ──► decimal. */
d2b: parse arg ?; return right( x2b( d2x(?) ), 64, 0) /*convert dec ──► 64 bit bin.*/
commas: parse arg _; do ?=length(_)-3 to 1 by -3; _= insert(',', _, ?); end; return _
/*──────────────────────────────────────────────────────────────────────────────────────*/
next: procedure expose state const o.; x= d2b(state) /*convert STATE to binary. */
x = xor(x, left( o.12 || x, 64) ) /*shift right 12 bits and XOR*/
x = xor(x, right( x || o.25, 64) ) /* " left 25 " " " */
x = xor(x, left( o.27 || x, 64) ) /* " right 27 " " " */
state= b2d(x) /*set STATE to the latest X.*/
return b2d( left( d2b(state * const), 32) ) /*return a pseudo─random num.*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
xor: parse arg a, b; $= /*perform a bit─wise XOR. */
do !=1 for length(a); $= $ || (substr(a,!,1) && substr(b,!,1) )
end /*!*/; return $
output   when using the default inputs:
 n      pseudo─random number
═══  ══════════════════════════
 1:       3,540,625,527
 2:       2,750,739,987
 3:       4,037,983,143
 4:       1,993,361,440
 5:       3,809,424,708

 #    count of pseudo─random #
═══  ══════════════════════════
 0:          20,103
 1:          19,922
 2:          19,937
 3:          20,031
 4:          20,007

Wren[edit]

Translation of: Python
Library: Wren-big

As Wren doesn't have a 64-bit integer type, we use BigInt instead.

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])")
Output:
3540625527
2750739987
4037983143
1993361440
3809424708

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