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

From Rosetta Code
Pseudo-random numbers/Splitmix64 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.

Splitmix64 is the default pseudo-random number generator algorithm in Java and is included / available in many other languages. It uses a fairly simple algorithm that, though it is considered to be poor for cryptographic purposes, is very fast to calculate, and is "good enough" for many random number needs. It passes several fairly rigorous PRNG "fitness" tests that some more complex algorithms fail.

Splitmix64 is not recommended for demanding random number requirements, but is often used to calculate initial states for other more complex pseudo-random number generators.

The "standard" splitmix64 maintains one 64 bit state variable and returns 64 bits of random data with each call.

Basic pseudocode algorithm:

    uint64 state                                  /* The state can be seeded with any (upto) 64 bit integer value. */

    next_int() {
        state += 0x9e3779b97f4a7c15               /* increment the state variable */
        uint64 z = state                          /* copy the state to a working variable */
        z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9  /* xor the variable with the variable right bit shifted 30 then multiply by a constant */
        z = (z ^ (z >> 27)) * 0x94d049bb133111eb  /* xor the variable with the variable right bit shifted 27 then multiply by a constant */
        return z ^ (z >> 31)                      /* return the variable xored with itself right bit shifted 31 */
    }

    next_float() {
        return next_int() / (1 << 64)             /* divide by 2^64 to return a value between 0 and 1 */
    }

The returned value should hold 64 bits of numeric data. If your language does not support unsigned 64 bit integers directly you may need to apply appropriate bitmasks during bitwise operations.

In keeping with the general layout of several recent pseudo-random number tasks:


Task
  • Write a class or set of functions that generates pseudo-random numbers using splitmix64.
  • Show the first five integers generated using the seed 1234567.
    6457827717110365317
    3203168211198807973
    9817491932198370423
    4593380528125082431
   16408922859458223821
  
  • Show that for an initial seed of 987654321, the counts of 100_000 repetitions of floor next_float() * 5 is as follows:
   0: 20027, 1: 19892, 2: 20073, 3: 19978, 4: 20030  
  • Show your output here, on this page.


See also


Related tasks


ALGOL 68[edit]

Works with: ALGOL 68G version Any Tested with release 2.8.3.win32
BEGIN # generate some pseudo random numbers using Splitmix64 #
# note that although LONG INT is 64 bits in Algol 68G, LONG BITS is longer than 64 bits #
LONG BITS mask 64 = LONG 16rffffffffffffffff;
LONG BITS state := 16r1234567;
LONG INT one shl 64 = ABS ( LONG 16r1 SHL 64 );
# 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 mask 64;
# add a LONG BITS value to a LONG BITS #
OP +:= = ( REF LONG BITS r, LONG BITS v )REF LONG BITS:
r := SHORTEN ( BIN ( LENG ABS r + LENG ABS v ) AND mask 64 );
# multiplies a LONG BITS value by a LONG BITS value #
OP *:= = ( REF LONG BITS r, LONG BITS v )REF LONG BITS:
r := SHORTEN ( BIN ( ABS LENG r * LENG ABS v ) AND mask 64 );
# gets the next pseudo random integer #
PROC next int = LONG INT:
BEGIN
state +:= LONG 16r9e3779b97f4a7c15;
LONG BITS z := state;
z XORAB ( z SHR 30 );
z *:= LONG 16rbf58476d1ce4e5b9;
z XORAB ( z SHR 27 );
z *:= LONG 16r94d049bb133111eb;
z XORAB ( z SHR 31 );
ABS z
END # next int # ;
# gets the next pseudo random real #
PROC next float = LONG REAL: next int / one shl 64;
BEGIN # task test cases #
seed( 1234567 );
print( ( whole( next int, 0 ), newline ) ); # 6457827717110365317 #
print( ( whole( next int, 0 ), newline ) ); # 3203168211198807973 #
print( ( whole( next int, 0 ), newline ) ); # 9817491932198370423 #
print( ( whole( next int, 0 ), newline ) ); # 4593380528125082431 #
print( ( whole( next int, 0 ), newline ) ); # 16408922859458223821 #
# 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:
6457827717110365317
3203168211198807973
9817491932198370423
4593380528125082431
16408922859458223821
 0:  20027 1:  19892 2:  20073 3:  19978 4:  20030

C[edit]

Code copied from the reference C implementation used by Java, and using GNU GCC v7.1.1.

/*  Written in 2015 by Sebastiano Vigna ([email protected])
 
To the extent possible under law, the author has dedicated all copyright
and related and neighboring rights to this software to the public domain
worldwide. This software is distributed without any warranty.
 
See <http://creativecommons.org/publicdomain/zero/1.0/>. */

 
#include <stdint.h>
#include <stdio.h>
#include <math.h>
 
/* This is a fixed-increment version of Java 8's SplittableRandom generator
See http://dx.doi.org/10.1145/2714064.2660195 and
http://docs.oracle.com/javase/8/docs/api/java/util/SplittableRandom.html
 
It is a very fast generator passing BigCrush, and it can be useful if
for some reason you absolutely want 64 bits of state. */

 
static uint64_t x; /* The state can be seeded with any value. */
 
uint64_t next() {
uint64_t z = (x += 0x9e3779b97f4a7c15);
z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
return z ^ (z >> 31);
}
 
double next_float() {
return next() / pow(2.0, 64);
}
 
int main() {
int i, j;
x = 1234567;
for(i = 0; i < 5; ++i)
printf("%llu\n", next()); /* needed to use %lu verb for GCC 7.5.0-3 */
x = 987654321;
int vec5[5] = {0, 0, 0, 0, 0};
for(i = 0; i < 100000; ++i) {
j = next_float() * 5.0;
vec5[j] += 1;
}
for(i = 0; i < 5; ++i)
printf("%d: %d ", i, vec5[i]);
}
 
Output:
6457827717110365317
3203168211198807973
9817491932198370423
4593380528125082431
16408922859458223821
0: 20027  1: 19892  2: 20073  3: 19978  4: 20030  

Factor[edit]

USING: io kernel math math.bitwise math.functions
math.statistics namespaces prettyprint sequences ;
 
SYMBOL: state
 
: seed ( n -- ) 64 bits state set ;
 
: next-int ( -- n )
0x9e3779b97f4a7c15 state [ + 64 bits ] change
state get -30 0xbf58476d1ce4e5b9 -27 0x94d049bb133111eb -31 1
[ [ dupd shift bitxor ] dip * 64 bits ] [email protected] ;
 
: next-float ( -- x ) next-int 64 2^ /f ;
 
! Test next-int
"Seed: 1234567; first five integer values" print
1234567 seed 5 [ next-int . ] times nl
 
! Test next-float
"Seed: 987654321; first 100,000 float values histogram" print
987654321 seed 100,000 [ next-float 5 * >integer ] replicate
histogram .
Output:
Seed: 1234567; first five integer values
6457827717110365317
3203168211198807973
9817491932198370423
4593380528125082431
16408922859458223821

Seed: 987654321; first 100,000 float values histogram
H{ { 0 20027 } { 1 19892 } { 2 20073 } { 3 19978 } { 4 20030 } }

Go[edit]

package main
 
import (
"fmt"
"math"
)
 
type Splitmix64 struct{ state uint64 }
 
func Splitmix64New(state uint64) *Splitmix64 { return &Splitmix64{state} }
 
func (sm64 *Splitmix64) nextInt() uint64 {
sm64.state += 0x9e3779b97f4a7c15
z := sm64.state
z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9
z = (z ^ (z >> 27)) * 0x94d049bb133111eb
return z ^ (z >> 31)
}
 
func (sm64 *Splitmix64) nextFloat() float64 {
return float64(sm64.nextInt()) / (1 << 64)
}
 
func main() {
randomGen := Splitmix64New(1234567)
for i := 0; i < 5; i++ {
fmt.Println(randomGen.nextInt())
}
 
var counts [5]int
randomGen = Splitmix64New(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:
6457827717110365317
3203168211198807973
9817491932198370423
4593380528125082431
16408922859458223821

The counts for 100,000 repetitions are:
  0 : 20027
  1 : 19892
  2 : 20073
  3 : 19978
  4 : 20030

Phix[edit]

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

include mpfr.e
mpz state = mpz_init(),
shift = mpz_init("0x9e3779b97f4a7c15"),
mult1 = mpz_init("0xbf58476d1ce4e5b9"),
mult2 = mpz_init("0x94d049bb133111eb"),
b64 = mpz_init("0x10000000000000000"), -- (truncate to 64 bits)
tmp = mpz_init(),
z = mpz_init()
 
procedure seed(integer num)
mpz_set_si(state,num)
end procedure
 
procedure next_int()
mpz_add(state, state, shift) -- state += shift
mpz_fdiv_r(state, state, b64) -- state := remainder(z,b64)
mpz_set(z, state) -- z := state
mpz_tdiv_q_2exp(tmp, z, 30) -- tmp := trunc(z/2^30)
mpz_xor(z, z, tmp) -- z := xor_bits(z,tmp)
mpz_mul(z, z, mult1) -- z *= mult1
mpz_fdiv_r(z, z, b64) -- z := remainder(z,b64)
mpz_tdiv_q_2exp(tmp, z, 27) -- tmp := trunc(z/2^27)
mpz_xor(z, z, tmp) -- z := xor_bits(z,tmp)
mpz_mul(z, z, mult2) -- z *= mult2
mpz_fdiv_r(z, z, b64) -- z := remainder(z,b64)
mpz_tdiv_q_2exp(tmp, z, 31) -- tmp := trunc(z/2^31)
mpz_xor(z, z, tmp) -- z := xor_bits(z,tmp)
-- (result left in z)
end procedure
 
function next_float()
next_int()
mpfr f = mpfr_init_set_z(z)
mpfr_div_z(f, f, b64)
return mpfr_get_d(f)
end function
 
seed(1234567)
for i=1 to 5 do
next_int()
printf(1,"%s\n",mpz_get_str(z))
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:
6457827717110365317
3203168211198807973
9817491932198370423
4593380528125082431
16408922859458223821
{20027,19892,20073,19978,20030}

Python[edit]

MASK64 = (1 << 64) - 1
C1 = 0x9e3779b97f4a7c15
C2 = 0xbf58476d1ce4e5b9
C3 = 0x94d049bb133111eb
 
 
 
class Splitmix64():
 
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**64"
z = self.state = (self.state + C1) & MASK64
z = ((z ^ (z >> 30)) * C2) & MASK64
z = ((z ^ (z >> 27)) * C3) & MASK64
answer = (z ^ (z >> 31)) & MASK64
 
return answer
 
def next_float(self):
"return random float between 0 and 1"
return self.next_int() / (1 << 64)
 
 
if __name__ == '__main__':
random_gen = Splitmix64()
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:
6457827717110365317
3203168211198807973
9817491932198370423
4593380528125082431
16408922859458223821
{0: 20027, 1: 19892, 2: 20073, 3: 19978, 4: 20030}

Raku[edit]

Works with: Rakudo version 2020.07
class splitmix64 {
has $!state;
 
submethod BUILD ( Int :$seed where * >= 0 = 1 ) { $!state = $seed }
 
method next-int {
my $next = $!state = ($!state + 0x9e3779b97f4a7c15) +& (2⁶⁴ - 1);
$next = ($next +^ ($next +> 30)) * 0xbf58476d1ce4e5b9 +& (2⁶⁴ - 1);
$next = ($next +^ ($next +> 27)) * 0x94d049bb133111eb +& (2⁶⁴ - 1);
($next +^ ($next +> 31)) +& (2⁶⁴ - 1);
}
 
method next-rat { self.next-int / 2⁶⁴ }
}
 
# Test next-int
say 'Seed: 1234567; first five Int values';
my $rng = splitmix64.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 = splitmix64.new( :seed(987654321) );
say ( ($rng.next-rat * 5).floor xx 100_000 ).Bag;
Output:
Seed: 1234567; first five Int values
6457827717110365317
3203168211198807973
9817491932198370423
4593380528125082431
16408922859458223821

Seed: 987654321; first 1e5 Rat values histogram
Bag(0(20027) 1(19892) 2(20073) 3(19978) 4(20030))

REXX[edit]

/*REXX program  generates   pseudo─random numbers   using the  split mix 64 bit  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.1= x2d( 9e3779b97f4a7c15 ) /*initialize 1st constant to be used. */
const.2= x2d('bf58476d1ce4e5b9') /* " 2nd " " " " */
const.3= x2d( 94d049bb133111eb ) /* " 3rd " " " " */
o.30= copies(0, 30) /*construct 30 bits of zeros. */
o.27= copies(0, 27) /* " 27 " " " */
o.31= copies(0, 31) /* " 31 " " " */
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()), 27) /*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**64 /*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(@.#), 15) /*show count of a random num.*/
end /*#*/
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg _; do ?=length(_)-3 to 1 by -3; _= insert(',', _, ?); end; return _
b2d: parse arg ?; return x2d( b2x(?) ) /*convert bin──►decimal. */
d2b: parse arg ?; return right( x2b( d2x(?) ), 64, 0) /*convert dec──►64 bit bin.*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
next: procedure expose state const. o.
state= state + const.1  ; z= d2b(state) /*add const1──►STATE; conv.*/
z= xor(z, left(o.30 || z, 64)); z= d2b(b2d(z)*const.2) /*shiftR 30 bits & XOR; " */
z= xor(z, left(o.27 || z, 64)); z= d2b(b2d(z)*const.3) /* " 27 " " " " */
z= xor(z, left(o.31 || z, 64)); return b2d(z) /* " 31 " " " " */
/*──────────────────────────────────────────────────────────────────────────────────────*/
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:    6,457,827,717,110,365,317
 2:    3,203,168,211,198,807,973
 3:    9,817,491,932,198,370,423
 4:    4,593,380,528,125,082,431
 5:   16,408,922,859,458,223,821

 #     count of pseudo─random #
═══  ════════════════════════════
 0:           20,027
 1:           19,892
 2:           20,073
 3:           19,978
 4:           20,030

Wren[edit]

Library: Wren-big

No 64 bit integers so we use BigInt with a mask.

import "/big" for BigInt
 
var Const1 = BigInt.fromBaseString("9e3779b97f4a7c15", 16)
var Const2 = BigInt.fromBaseString("bf58476d1ce4e5b9", 16)
var Const3 = BigInt.fromBaseString("94d049bb133111eb", 16)
var Mask64 = (BigInt.one << 64) - BigInt.one
 
class Splitmix64 {
construct new(state) {
_state = state
}
 
nextInt {
_state = (_state + Const1) & Mask64
var z = _state
z = ((z ^ (z >> 30)) * Const2) & Mask64
z = ((z ^ (z >> 27)) * Const3) & Mask64
return (z ^ (z >> 31)) & Mask64
}
 
nextFloat { nextInt.toNum / 2.pow(64) }
}
 
var randomGen = Splitmix64.new(BigInt.new(1234567))
for (i in 0..4) System.print(randomGen.nextInt)
 
var counts = List.filled(5, 0)
randomGen = Splitmix64.new(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:
6457827717110365317
3203168211198807973
9817491932198370423
4593380528125082431
16408922859458223821

The counts for 100,000 repetitions are:
  0 : 20027
  1 : 19892
  2 : 20073
  3 : 19978
  4 : 20030