Pseudo-random numbers/Splitmix64: Difference between revisions

m
→‎jq: 13 minutes u+s
No edit summary
m (→‎jq: 13 minutes u+s)
(4 intermediate revisions by 4 users not shown)
Line 732:
}
exec test(0x1000:8)</syntaxhighlight>
 
=={{header|jq}}==
'''Works with gojq, the Go implementation of jq'''
 
'''Adapted from [[#Wren|Wren]]'''
 
This entry assumes sufficiently precise integer arithmetic, e.g. as provided by gojq,
which supports infinite-precision integer arithmetic.
Unfortunately, though, gojq does not support efficient bitwise operations, and using
the following to generate 100,000 random numbers takes about 13 minutes on a 3GHz machine.
 
The main significance of this entry is thus to increase confidence in the correctness of the
definitions as well as in the Go implementation of jq.
 
In the following, a 'bitarray' is a 0/1 array, and when integers are represented by bitarrays, the first bit is the least-significant one. Unless otherwise indicated, the bitwise operations work on bitarrays.
 
'''Generic utilities'''
<syntaxhighlight lang="jq">
# Input: a string in base $b (2 to 35 inclusive)
# Output: a JSON number, being the decimal value corresponding to the input.
def frombase($b):
def decimalValue:
if 48 <= . and . <= 57 then . - 48
elif 65 <= . and . <= 90 then . - 55 # (10+.-65)
elif 97 <= . and . <= 122 then . - 87 # (10+.-97)
else "decimalValue" | error
end;
reduce (explode|reverse[]|decimalValue) as $x ({p:1};
.value += (.p * $x)
| .p *= $b)
| .value ;
 
# To take advantage of gojq's arbitrary-precision integer arithmetic:
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in);
 
# If the input and $j are integers, then the result will be an integer.
def div($j):
(. - (. % j)) / $j;
 
# Convert an integer to a bitarray, least significant bit first
def bitwise:
recurse( if . >= 2 then div(2) else empty end) | . % 2;
 
# Essentially the inverse of bitwise,
# i.e. interpret an array of 0s and 1s (with least-significant-bit first ) as a decimal
def to_int:
. as $in
# state: [sum, power]
| reduce .[] as $i ([0, 1]; .[1] as $p | [.[0] + $p * $i, ($p * 2)])
| .[0];
 
# $x and $y and output are bitarrays
def xor($x;$y):
def lxor(a;b):
if (a==1 or b==1) and ((a==1 and b==1)|not) then 1
elif a == null then b
elif b == null then a
else 0
end;
if $x == [0] then $y
elif $y == [0] then $x
else
[ range(0; [($x|length), ($y|length)] | max) as $i
| lxor($x[$i]; $y[$i]) ]
end ;
 
# $x and $y and output are bitarrays
def xand($x;$y):
def lxand(a;b):
(a==1 and b==1) | 1 // 0;
if $x == [0] or $y == [0] then [0]
else
[range(0; [($x|length), ($y|length)] | min) as $i
| lxand($x[$i]; $y[$i]) ]
end ;
 
# shift right
def right($n): .[$n:];
 
def mask64: .[:64];
 
# input and output: a bitarray
def mult($int):
($int * to_int) | [bitwise];
 
def plus($int):
($int + to_int) | [bitwise];
 
def tabulate(stream):
reduce stream as $i ([]; .[$i] += 1)
| range(0;length) as $i
| " \($i) : \(.[$i] // 0)" ;
</syntaxhighlight>
'''Splitmix64'''
<syntaxhighlight lang="jq">
# input: a bitarray
def nextInt:
def Const1: "9e3779b97f4a7c15" | frombase(16) ;
def Const2: "bf58476d1ce4e5b9" | frombase(16) ;
def Const3: "94d049bb133111eb" | frombase(16) ;
 
(plus(Const1) | mask64)
| . as $state
| xor(.; right(30)) | mult(Const2) | mask64
| xor(.; right(27)) | mult(Const3) | mask64
| xor(.; right(31)) | mask64
| ., ($state|nextInt) ;
 
def randomInt64: [bitwise] | nextInt | to_int;
 
def randomReal:
pow(2;64) as $d
| [bitwise] | nextInt | to_int / $d;
 
### The tasks
(limit(5; 1234567 | randomInt64)),
 
"\nThe counts for 100,000 repetitions are:",
tabulate( limit(100; 987654321 | randomReal * 5 | floor) )
 
</syntaxhighlight>
{{output}}
<pre>
6457827717110365317
3203168211198807973
9817491932198370423
4593380528125082431
16408922859458223821
 
The counts for 100,000 repetitions are:
0 : 20027
1 : 19892
2 : 20073
3 : 19978
4 : 20030
</pre>
 
=={{header|Julia}}==
Line 1,049 ⟶ 1,185:
=={{header|Phix}}==
As per [[Pseudo-random_numbers/PCG32#Phix]], resorting to mpfr/gmp
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
include mpfr.e
<span style="color: #008080;">include</span> <span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span>
mpz state = mpz_init(),
<span style="color: #004080;">mpz</span> <span style="color: #000000;">state</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(),</span>
shift = mpz_init("0x9e3779b97f4a7c15"),
<span style="color: #000000;">shift</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"0x9e3779b97f4a7c15"</span><span style="color: #0000FF;">),</span>
mult1 = mpz_init("0xbf58476d1ce4e5b9"),
<span style="color: #000000;">mult1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"0xbf58476d1ce4e5b9"</span><span style="color: #0000FF;">),</span>
mult2 = mpz_init("0x94d049bb133111eb"),
<span style="color: #000000;">mult2</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"0x94d049bb133111eb"</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">b64</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"0x10000000000000000"</span><span style="color: #0000FF;">),</span> <span style="color: #000080;font-style:italic;">-- (truncate to 64 bits)</span>
tmp = mpz_init(),
<span style="color: #000000;">tmp</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(),</span>
z = mpz_init()
<span style="color: #000000;">z</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">()</span>
procedure seed(integer num)
<span style="color: #008080;">procedure</span> <span style="color: #000000;">seed</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">num</span><span style="color: #0000FF;">)</span>
mpz_set_si(state,num)
<span style="color: #7060A8;">mpz_set_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">state</span><span style="color: #0000FF;">,</span><span style="color: #000000;">num</span><span style="color: #0000FF;">)</span>
end procedure
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
procedure next_int()
<span style="color: #008080;">procedure</span> <span style="color: #000000;">next_int</span><span style="color: #0000FF;">()</span>
mpz_add(state, state, shift) -- state += shift
<span style="color: #7060A8;">mpz_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">state</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">state</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">shift</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- state += shift</span>
mpz_fdiv_r(state, state, b64) -- state := remainder(z,b64)
<span style="color: #7060A8;">mpz_fdiv_r</span><span style="color: #0000FF;">(</span><span style="color: #000000;">state</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">state</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b64</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- state := remainder(z,b64) </span>
mpz_set(z, state) -- z := state
<span style="color: #7060A8;">mpz_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">state</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- z := state</span>
mpz_tdiv_q_2exp(tmp, z, 30) -- tmp := trunc(z/2^30)
<span style="color: #7060A8;">mpz_tdiv_q_2exp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tmp</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">30</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- tmp := trunc(z/2^30)</span>
mpz_xor(z, z, tmp) -- z := xor_bits(z,tmp)
<span style="color: #7060A8;">mpz_xor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">tmp</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- z := xor_bits(z,tmp)</span>
mpz_mul(z, z, mult1) -- z *= mult1
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">mult1</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- z *= mult1</span>
mpz_fdiv_r(z, z, b64) -- z := remainder(z,b64)
<span style="color: #7060A8;">mpz_fdiv_r</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b64</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- z := remainder(z,b64) </span>
mpz_tdiv_q_2exp(tmp, z, 27) -- tmp := trunc(z/2^27)
<span style="color: #7060A8;">mpz_tdiv_q_2exp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tmp</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">27</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- tmp := trunc(z/2^27)</span>
mpz_xor(z, z, tmp) -- z := xor_bits(z,tmp)
<span style="color: #7060A8;">mpz_xor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">tmp</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- z := xor_bits(z,tmp)</span>
mpz_mul(z, z, mult2) -- z *= mult2
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">mult2</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- z *= mult2</span>
mpz_fdiv_r(z, z, b64) -- z := remainder(z,b64)
<span style="color: #7060A8;">mpz_fdiv_r</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b64</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- z := remainder(z,b64) </span>
mpz_tdiv_q_2exp(tmp, z, 31) -- tmp := trunc(z/2^31)
<span style="color: #7060A8;">mpz_tdiv_q_2exp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tmp</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">31</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- tmp := trunc(z/2^31)</span>
mpz_xor(z, z, tmp) -- z := xor_bits(z,tmp)
<span style="color: #7060A8;">mpz_xor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">tmp</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- z := xor_bits(z,tmp)
-- (result left in z)</span>
end procedure
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
function next_float()
<span style="color: #008080;">function</span> <span style="color: #000000;">next_float</span><span style="color: #0000FF;">()</span>
next_int()
<span style="color: #000000;">next_int</span><span style="color: #0000FF;">()</span>
mpfr f = mpfr_init_set_z(z)
<span style="color: #004080;">mpfr</span> <span style="color: #000000;">f</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init_set_z</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)</span>
mpfr_div_z(f, f, b64)
<span style="color: #7060A8;">mpfr_div_z</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b64</span><span style="color: #0000FF;">)</span>
return mpfr_get_d(f)
<span style="color: #008080;">return</span> <span style="color: #7060A8;">mpfr_get_d</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">)</span>
end function
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
seed(1234567)
<span style="color: #000000;">seed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1234567</span><span style="color: #0000FF;">)</span>
for i=1 to 5 do
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">5</span> <span style="color: #008080;">do</span>
next_int()
<span style="color: #000000;">next_int</span><span style="color: #0000FF;">()</span>
printf(1,"%s\n",mpz_get_str(z))
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%s\n"</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">))</span>
end for
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
seed(987654321)
<span style="color: #000000;">seed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">987654321</span><span style="color: #0000FF;">)</span>
sequence r = repeat(0,5)
<span style="color: #004080;">sequence</span> <span style="color: #000000;">r</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">5</span><span style="color: #0000FF;">)</span>
for i=1 to 100000 do
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">100000</span> <span style="color: #008080;">do</span>
integer rdx = floor(next_float()*5)+1
<span style="color: #004080;">integer</span> <span style="color: #000000;">rdx</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">next_float</span><span style="color: #0000FF;">()*</span><span style="color: #000000;">5</span><span style="color: #0000FF;">)+</span><span style="color: #000000;">1</span>
r[rdx] += 1
<span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">rdx</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
end for
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
?r
<span style="color: #0000FF;">?</span><span style="color: #000000;">r</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
Line 1,363 ⟶ 1,499:
3: 19,978
4: 20,030
</pre>
 
=={{header|RPL}}==
RPL offers hundreds of instructions and commands but, strangely, bitwise operations are limited to the set defined for the HP-16C RPN calculator in 1982.
So, to shift an integer n times, we can either stay in low gear, shifting one bit at a time:
≪ 1 SWAP '''START''' SR '''NEXT'''
≫ '<span style="color:blue">SRn</span>' STO
or play with the gearbox:
8 MOD LAST / IP ROT <span style="color:grey">@ q, r = divmod(n,8)</span>
'''IF''' SWAP THEN 1 LAST '''START''' SRB '''NEXT END''' <span style="color:grey">@ q, r = divmod(n,8)</span>
'''IF''' SWAP '''THEN''' 1 LAST '''START''' SR '''NEXT END''' <span style="color:grey">@ shift & bit right, r times</span>
≫ '<span style="color:blue">SRn</span>' STO
The first version is less efficient in terms of response time, but more idiomatic because it takes much less time to type on a calculator keyboard.
{{works with|HP|28}}
≪ # 9E3779B97F4A7C15h 'STATE' STO+
<span style="color:green">STATE</span> DUP 30 SRn XOR # BF58476D1CE4E5B9h *
DUP 27 <span style="color:blue">SRn</span> XOR # 94D049BB133111EBh *
DUP 31 <span style="color:blue">SRn</span> XOR
≫ '<span style="color:blue">NEXTINT</span>' STO
≪ <span style="color:blue">NEXTINT</span> B→R 2 64 ^ /
≫ '<span style="color:blue">NEXTFLOAT</span>' STO
≪ # 1234567d '<span style="color:green">STATE</span>' STO
1 5 '''START''' <span style="color:blue">NEXTINT</span> '''NEXT'''
≫ '<span style="color:blue">TASK1</span>' STO
≪ { 5 } O CON
1 10000 '''START'''
<span style="color:blue">NEXTFLOAT</span> 5 * CEIL DUP2 GET 1 + PUT '''NEXT'''
≫ '<span style="color:blue">TASK2</span>' STO
 
The last task requirement takes too much time to be run on a calculator. On a emulator, we had to split it into 10 runs to not wake up the watchdog timer.
<span style="color:blue">TASK1</span>
# 987654321d '<span style="color:green">STATE</span>' STO
{ 5 } O CON
<span style="color:blue">TASK2</span> <span style="color:blue">TASK2</span> .. <span style="color:blue">TASK2</span> <span style="color:grey">@ 10 times</span>
{{out}}
<pre>
6: # 6457827717110365317d
5: # 3203168211198807973d
4: # 9817491932198370423d
3: # 4593380528125082431d
2: # 16408922859458223821d
1: [ 20027 19892 20073 19978 20030 ]
</pre>
 
Line 1,397 ⟶ 1,579:
{0=>20027, 1=>19892, 2=>20073, 3=>19978, 4=>20030}
</pre>
 
=={{header|Sidef}}==
{{trans|Perl}}
Line 1,445 ⟶ 1,628:
{{libheader|Wren-big}}
No 64 bit integers so we use BigInt with a mask.
<syntaxhighlight lang="ecmascriptwren">import "./big" for BigInt
 
var Const1 = BigInt.fromBaseString("9e3779b97f4a7c15", 16)
2,442

edits