Pseudo-random numbers/Xorshift star: Difference between revisions
m (→{{header|Julia}}: no masking for 64 bit unsigneds) |
(Added Algol 68) |
||
Line 82: | Line 82: | ||
* Show your output here, on this page. |
* Show your output here, on this page. |
||
=={{header|ALGOL 68}}== |
|||
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}} |
|||
<lang algol68>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</lang> |
|||
{{out}} |
|||
<pre> |
|||
3540625527 |
|||
2750739987 |
|||
4037983143 |
|||
1993361440 |
|||
3809424708 |
|||
0: 20103 1: 19922 2: 19937 3: 20031 4: 20007 |
|||
</pre> |
|||
=={{header|F_Sharp|F#}}== |
=={{header|F_Sharp|F#}}== |
||
Line 110: | Line 161: | ||
(4, 20007)(2, 19937)(3, 20031)(0, 20103)(1, 19922) |
(4, 20007)(2, 19937)(3, 20031)(0, 20103)(1, 19922) |
||
</pre> |
</pre> |
||
=={{header|Factor}}== |
=={{header|Factor}}== |
||
<lang factor>USING: accessors kernel literals math math.statistics |
<lang factor>USING: accessors kernel literals math math.statistics |
Revision as of 09:21, 15 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
- 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
<lang algol68>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</lang>
- Output:
3540625527 2750739987 4037983143 1993361440 3809424708 0: 20103 1: 19922 2: 19937 3: 20031 4: 20007
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
<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
<lang julia>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()
</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
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³² }
}
- 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
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