Pseudo-random numbers/Combined recursive generator MRG32k3a: Difference between revisions
(22 intermediate revisions by 11 users not shown) | |||
Line 73: | Line 73: | ||
{{trans|Python}} |
{{trans|Python}} |
||
< |
<syntaxhighlight lang="11l">V a1 = [Int64(0), 1403580, -810728] |
||
V m1 = Int64(2) ^ 32 - 209 |
V m1 = Int64(2) ^ 32 - 209 |
||
V a2 = [Int64(527612), 0, -1370589] |
V a2 = [Int64(527612), 0, -1370589] |
||
Line 112: | Line 112: | ||
L 100'000 |
L 100'000 |
||
hist[Int(random_gen.next_float() * 5)]++ |
hist[Int(random_gen.next_float() * 5)]++ |
||
print(hist)</ |
print(hist)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 125: | Line 125: | ||
=={{header|Ada}}== |
=={{header|Ada}}== |
||
< |
<syntaxhighlight lang="ada">package MRG32KA is |
||
type I64 is range -2**63..2**63 - 1; |
type I64 is range -2**63..2**63 - 1; |
||
m1 : constant I64 := 2**32 - 209; |
m1 : constant I64 := 2**32 - 209; |
||
Line 136: | Line 136: | ||
function Next_Float return Long_Float; |
function Next_Float return Long_Float; |
||
end MRG32KA; |
end MRG32KA; |
||
</syntaxhighlight> |
|||
</lang> |
|||
<syntaxhighlight lang="ada"> |
|||
<lang Ada> |
|||
package body MRG32KA is |
package body MRG32KA is |
||
Line 192: | Line 192: | ||
end MRG32KA; |
end MRG32KA; |
||
</syntaxhighlight> |
|||
</lang> |
|||
< |
<syntaxhighlight lang="ada">with Ada.Text_IO; use Ada.Text_IO; |
||
with mrg32ka; use mrg32ka; |
with mrg32ka; use mrg32ka; |
||
Line 219: | Line 219: | ||
end Main; |
end Main; |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{output}} |
{{output}} |
||
<pre> |
<pre> |
||
Line 232: | Line 232: | ||
=={{header|C}}== |
=={{header|C}}== |
||
< |
<syntaxhighlight lang="c">#include <math.h> |
||
#include <stdio.h> |
#include <stdio.h> |
||
#include <stdint.h> |
#include <stdint.h> |
||
Line 317: | Line 317: | ||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
|||
<pre>1459213977 |
|||
2827710106 |
|||
4245671317 |
|||
3877608661 |
|||
2595287583 |
|||
0: 20002 |
|||
1: 20060 |
|||
2: 19948 |
|||
3: 20059 |
|||
4: 19931</pre> |
|||
=={{header|C++}}== |
|||
{{trans|C}} |
|||
<syntaxhighlight lang="cpp">#include <array> |
|||
#include <iostream> |
|||
int64_t mod(int64_t x, int64_t y) { |
|||
int64_t m = x % y; |
|||
if (m < 0) { |
|||
if (y < 0) { |
|||
return m - y; |
|||
} else { |
|||
return m + y; |
|||
} |
|||
} |
|||
return m; |
|||
} |
|||
class RNG { |
|||
private: |
|||
// First generator |
|||
const std::array<int64_t, 3> a1{ 0, 1403580, -810728 }; |
|||
const int64_t m1 = (1LL << 32) - 209; |
|||
std::array<int64_t, 3> x1; |
|||
// Second generator |
|||
const std::array<int64_t, 3> a2{ 527612, 0, -1370589 }; |
|||
const int64_t m2 = (1LL << 32) - 22853; |
|||
std::array<int64_t, 3> x2; |
|||
// other |
|||
const int64_t d = (1LL << 32) - 209 + 1; // m1 + 1 |
|||
public: |
|||
void seed(int64_t state) { |
|||
x1 = { state, 0, 0 }; |
|||
x2 = { state, 0, 0 }; |
|||
} |
|||
int64_t next_int() { |
|||
int64_t x1i = mod((a1[0] * x1[0] + a1[1] * x1[1] + a1[2] * x1[2]), m1); |
|||
int64_t x2i = mod((a2[0] * x2[0] + a2[1] * x2[1] + a2[2] * x2[2]), m2); |
|||
int64_t z = mod(x1i - x2i, m1); |
|||
// keep last three values of the first generator |
|||
x1 = { x1i, x1[0], x1[1] }; |
|||
// keep last three values of the second generator |
|||
x2 = { x2i, x2[0], x2[1] }; |
|||
return z + 1; |
|||
} |
|||
double next_float() { |
|||
return static_cast<double>(next_int()) / d; |
|||
} |
|||
}; |
|||
int main() { |
|||
RNG rng; |
|||
rng.seed(1234567); |
|||
std::cout << rng.next_int() << '\n'; |
|||
std::cout << rng.next_int() << '\n'; |
|||
std::cout << rng.next_int() << '\n'; |
|||
std::cout << rng.next_int() << '\n'; |
|||
std::cout << rng.next_int() << '\n'; |
|||
std::cout << '\n'; |
|||
std::array<int, 5> counts{ 0, 0, 0, 0, 0 }; |
|||
rng.seed(987654321); |
|||
for (size_t i = 0; i < 100000; i++) { |
|||
auto value = floor(rng.next_float() * 5.0); |
|||
counts[value]++; |
|||
} |
|||
for (size_t i = 0; i < counts.size(); i++) { |
|||
std::cout << i << ": " << counts[i] << '\n'; |
|||
} |
|||
return 0; |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>1459213977 |
|||
2827710106 |
|||
4245671317 |
|||
3877608661 |
|||
2595287583 |
|||
0: 20002 |
|||
1: 20060 |
|||
2: 19948 |
|||
3: 20059 |
|||
4: 19931</pre> |
|||
=={{header|D}}== |
|||
{{trans|C++}} |
|||
<syntaxhighlight lang="d">import std.math; |
|||
import std.stdio; |
|||
long mod(long x, long y) { |
|||
long m = x % y; |
|||
if (m < 0) { |
|||
if (y < 0) { |
|||
return m - y; |
|||
} else { |
|||
return m + y; |
|||
} |
|||
} |
|||
return m; |
|||
} |
|||
class RNG { |
|||
private: |
|||
// First generator |
|||
immutable(long []) a1 = [0, 1403580, -810728]; |
|||
immutable long m1 = (1L << 32) - 209; |
|||
long[3] x1; |
|||
// Second generator |
|||
immutable(long []) a2 = [527612, 0, -1370589]; |
|||
immutable long m2 = (1L << 32) - 22853; |
|||
long[3] x2; |
|||
// other |
|||
immutable long d = m1 + 1; |
|||
public: |
|||
void seed(long state) { |
|||
x1 = [state, 0, 0]; |
|||
x2 = [state, 0, 0]; |
|||
} |
|||
long next_int() { |
|||
long x1i = mod((a1[0] * x1[0] + a1[1] * x1[1] + a1[2] * x1[2]), m1); |
|||
long x2i = mod((a2[0] * x2[0] + a2[1] * x2[1] + a2[2] * x2[2]), m2); |
|||
long z = mod(x1i - x2i, m1); |
|||
// keep the last three values of the first generator |
|||
x1 = [x1i, x1[0], x1[1]]; |
|||
// keep the last three values of the second generator |
|||
x2 = [x2i, x2[0], x2[1]]; |
|||
return z + 1; |
|||
} |
|||
double next_float() { |
|||
return cast(double) next_int() / d; |
|||
} |
|||
} |
|||
void main() { |
|||
auto rng = new RNG(); |
|||
rng.seed(1234567); |
|||
writeln(rng.next_int); |
|||
writeln(rng.next_int); |
|||
writeln(rng.next_int); |
|||
writeln(rng.next_int); |
|||
writeln(rng.next_int); |
|||
writeln; |
|||
int[5] counts; |
|||
rng.seed(987654321); |
|||
foreach (i; 0 .. 100_000) { |
|||
auto value = cast(int) floor(rng.next_float * 5.0); |
|||
counts[value]++; |
|||
} |
|||
foreach (i,v; counts) { |
|||
writeln(i, ": ", v); |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
{{out}} |
||
<pre>1459213977 |
<pre>1459213977 |
||
Line 332: | Line 510: | ||
=={{header|Factor}}== |
=={{header|Factor}}== |
||
< |
<syntaxhighlight lang="factor">USING: arrays kernel math math.order math.statistics |
||
math.vectors prettyprint sequences ; |
math.vectors prettyprint sequences ; |
||
Line 357: | Line 535: | ||
987654321 seed 100,000 [ next-float 5 * >integer ] replicate |
987654321 seed 100,000 [ next-float 5 * >integer ] replicate |
||
2nip histogram .</ |
2nip histogram .</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 368: | Line 546: | ||
</pre> |
</pre> |
||
=={{header|Forth}}== |
|||
{{trans|uBasic/4tH}} |
|||
{{works with|4tH v3.64}} |
|||
<syntaxhighlight lang="forth">6 array (seed) \ holds the seed |
|||
6 array (gens) \ holds the generators |
|||
\ set up constants |
|||
0 (gens) 0 th ! \ 1st generator |
|||
1403580 (gens) 1 th ! |
|||
-810728 (gens) 2 th ! |
|||
527612 (gens) 3 th ! \ 2nd generator |
|||
0 (gens) 4 th ! |
|||
-1370589 (gens) 5 th ! |
|||
1 32 lshift 209 - value (m) \ 1st generator constant |
|||
1 32 lshift 22853 - value (n) \ 2nd generator constant |
|||
( n1 n2 -- n3) |
|||
: (mod) tuck mod tuck 0< if abs + ;then drop ; |
|||
: (generate) do (seed) i th @ (gens) i th @ * + loop swap (mod) ; |
|||
: (reseed) ?do (seed) i th ! loop ; ( n1 n2 n3 limit index --) |
|||
: randomize 6 0 do dup i 3 mod if >zero then (seed) i th ! loop drop ; |
|||
( n --) |
|||
: random ( -- n) |
|||
(m) 0 3 0 (generate) (n) 0 6 3 (generate) over over |
|||
(seed) 4 th @ (seed) 3 th @ rot 6 3 (reseed) |
|||
(seed) 1 th @ (seed) 0 th @ rot 3 0 (reseed) - (m) (mod) 1+ |
|||
; |
|||
include lib/fp1.4th \ simple floating point support |
|||
include lib/zenfloor.4th \ for FLOOR |
|||
5 array (count) \ setup an array of 5 elements |
|||
: test |
|||
1234567 randomize |
|||
random . cr random . cr random . cr |
|||
random . cr random . cr cr \ perform the first test |
|||
987654321 randomize (m) 1+ s>f \ set up denominator |
|||
100000 0 ?do \ do this 100,000 times |
|||
random s>f fover f/ 5 s>f f* floor f>s cells (count) + 1 swap +! |
|||
loop fdrop |
|||
\ show the results |
|||
5 0 ?do i . ." : " (count) i th ? cr loop |
|||
; |
|||
test</syntaxhighlight> |
|||
{{out}} |
|||
<pre>1459213977 |
|||
2827710106 |
|||
4245671317 |
|||
3877608661 |
|||
2595287583 |
|||
0 : 20002 |
|||
1 : 20060 |
|||
2 : 19948 |
|||
3 : 20059 |
|||
4 : 19931 |
|||
</pre> |
|||
=={{header|Go}}== |
=={{header|Go}}== |
||
{{trans|Python}} |
{{trans|Python}} |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 437: | Line 675: | ||
fmt.Printf(" %d : %d\n", i, counts[i]) |
fmt.Printf(" %d : %d\n", i, counts[i]) |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 448: | Line 686: | ||
The counts for 100,000 repetitions are: |
The counts for 100,000 repetitions are: |
||
0 : 20002 |
|||
1 : 20060 |
|||
2 : 19948 |
|||
3 : 20059 |
|||
4 : 19931 |
|||
</pre> |
|||
=={{header|Haskell}}== |
|||
<syntaxhighlight lang="haskell">import Data.List |
|||
randoms :: Int -> [Int] |
|||
randoms seed = unfoldr go ([seed,0,0],[seed,0,0]) |
|||
where |
|||
go (x1,x2) = |
|||
let x1i = sum (zipWith (*) x1 a1) `mod` m1 |
|||
x2i = sum (zipWith (*) x2 a2) `mod` m2 |
|||
in Just $ ((x1i - x2i) `mod` m1, (x1i:init x1, x2i:init x2)) |
|||
a1 = [0, 1403580, -810728] |
|||
m1 = 2^32 - 209 |
|||
a2 = [527612, 0, -1370589] |
|||
m2 = 2^32 - 22853 |
|||
randomsFloat = map ((/ (2^32 - 208)) . fromIntegral) . randoms</syntaxhighlight> |
|||
<pre>*Main> take 5 $ randoms 1234567 |
|||
[1459213976,2827710105,4245671316,3877608660,2595287582] |
|||
*Main> let hist = map length . group . sort |
|||
*Main> hist . take 100000 $ (floor . (*5)) <$> randomsFloat 987654321 |
|||
[20002,20060,19948,20059,19931]</pre> |
|||
=== As a RandomGen instanse === |
|||
<syntaxhighlight lang="haskell">import System.Random |
|||
newtype MRG32k3a = MRG32k3a ([Int],[Int]) |
|||
mkMRG32k3a s = MRG32k3a ([s,0,0],[s,0,0]) |
|||
instance RandomGen MRG32k3a where |
|||
next (MRG32k3a (x1,x2)) = |
|||
let x1i = sum (zipWith (*) x1 a1) `mod` m1 |
|||
x2i = sum (zipWith (*) x2 a2) `mod` m2 |
|||
in ((x1i - x2i) `mod` m1, MRG32k3a (x1i:init x1, x2i:init x2)) |
|||
where |
|||
a1 = [0, 1403580, -810728] |
|||
m1 = 2^32 - 209 |
|||
a2 = [527612, 0, -1370589] |
|||
m2 = 2^32 - 22853 |
|||
split _ = error "MRG32k3a is not splittable"</syntaxhighlight> |
|||
In this case the sequence or numbers differs from direct unfolding, due to internal uniform shuffling. |
|||
<pre>*Main> take 5 $ randoms (mkMRG32k3a 1234567) |
|||
[2827710105,3877608660,3642754129,1259674122,3002249941] |
|||
*Main> let hist = map length . group . sort |
|||
*Main> hist . take 100000 $ (floor . (*5)) <$> (randoms (mkMRG32k3a 987654321) :: [Float]) |
|||
[20015,19789,20024,20133,20039]</pre> |
|||
=={{header|Java}}== |
|||
{{trans|C++}} |
|||
<syntaxhighlight lang="java">public class App { |
|||
private static long mod(long x, long y) { |
|||
long m = x % y; |
|||
if (m < 0) { |
|||
if (y < 0) { |
|||
return m - y; |
|||
} else { |
|||
return m + y; |
|||
} |
|||
} |
|||
return m; |
|||
} |
|||
public static class RNG { |
|||
// first generator |
|||
private final long[] a1 = {0, 1403580, -810728}; |
|||
private static final long m1 = (1L << 32) - 209; |
|||
private long[] x1; |
|||
// second generator |
|||
private final long[] a2 = {527612, 0, -1370589}; |
|||
private static final long m2 = (1L << 32) - 22853; |
|||
private long[] x2; |
|||
// other |
|||
private static final long d = m1 + 1; |
|||
public void seed(long state) { |
|||
x1 = new long[]{state, 0, 0}; |
|||
x2 = new long[]{state, 0, 0}; |
|||
} |
|||
public long nextInt() { |
|||
long x1i = mod(a1[0] * x1[0] + a1[1] * x1[1] + a1[2] * x1[2], m1); |
|||
long x2i = mod(a2[0] * x2[0] + a2[1] * x2[1] + a2[2] * x2[2], m2); |
|||
long z = mod(x1i - x2i, m1); |
|||
// keep the last three values of the first generator |
|||
x1 = new long[]{x1i, x1[0], x1[1]}; |
|||
// keep the last three values of the second generator |
|||
x2 = new long[]{x2i, x2[0], x2[1]}; |
|||
return z + 1; |
|||
} |
|||
public double nextFloat() { |
|||
return 1.0 * nextInt() / d; |
|||
} |
|||
} |
|||
public static void main(String[] args) { |
|||
RNG rng = new RNG(); |
|||
rng.seed(1234567); |
|||
System.out.println(rng.nextInt()); |
|||
System.out.println(rng.nextInt()); |
|||
System.out.println(rng.nextInt()); |
|||
System.out.println(rng.nextInt()); |
|||
System.out.println(rng.nextInt()); |
|||
System.out.println(); |
|||
int[] counts = {0, 0, 0, 0, 0}; |
|||
rng.seed(987654321); |
|||
for (int i = 0; i < 100_000; i++) { |
|||
int value = (int) Math.floor(rng.nextFloat() * 5.0); |
|||
counts[value]++; |
|||
} |
|||
for (int i = 0; i < counts.length; i++) { |
|||
System.out.printf("%d: %d%n", i, counts[i]); |
|||
} |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>1459213977 |
|||
2827710106 |
|||
4245671317 |
|||
3877608661 |
|||
2595287583 |
|||
0: 20002 |
|||
1: 20060 |
|||
2: 19948 |
|||
3: 20059 |
|||
4: 19931</pre> |
|||
=={{header|jq}}== |
|||
'''Adapted from [[#Wren|Wren]]''' |
|||
'''Works with jq, the C implementation of jq''' |
|||
'''Works with gojq, the Go implementation of jq''' |
|||
'''Works with jaq, the Rust implementation of jq''' |
|||
<syntaxhighlight lang="jq"> |
|||
# constants |
|||
def A1: [0, 1403580, -810728]; |
|||
def M1: 4294967087; # pow(2;32) - 209 |
|||
def A2: [527612, 0, -1370589]; |
|||
def M2: 4294944443; # pow(2;32) - 22853 |
|||
def D: M1 + 1; |
|||
# Python-style modulus |
|||
def Mod($x; $y): |
|||
def abs: if . < 0 then - . else . end; |
|||
($x % $y) as $m |
|||
| if $m < 0 then $m + ($y|abs) else $m end; |
|||
def MRG32k3a: {x1: [0, 0, 0], x2: [0, 0, 0] }; |
|||
def seed($seedState): |
|||
if $seedState <= 0 or $seedState >= D |
|||
then "Argument of seed must be in the range (0, \(D))" | error |
|||
else {x1: [$seedState, 0, 0]} | .x2 = .x1 |
|||
end ; |
|||
# Input: {x1, x2} as for MRG32k31 |
|||
# Output: {x1, x2, nextInt} |
|||
def nextInt: |
|||
Mod(A1[0]*.x1[0] + A1[1]*.x1[1] + A1[2]*.x1[2]; M1) as $x1i |
|||
| Mod(A2[0]*.x2[0] + A2[1]*.x2[1] + A2[2]*.x2[2]; M2) as $x2i |
|||
| .x1 = [$x1i, .x1[0], .x1[1]] # keep last three |
|||
| .x2 = [$x2i, .x2[0], .x2[1]] # keep last three |
|||
| .nextInt = Mod($x1i - $x2i; M1) + 1 ; |
|||
def nextFloat: nextInt | .nextFloat = (.nextInt / D); |
|||
def task1: |
|||
foreach range(0; 5) as $i (seed(1234567); nextInt ) | .nextInt; |
|||
def task2($n): |
|||
seed(987654321) |
|||
| reduce range(0; $n) as $i (.counts = [range(0;5)|0]; |
|||
nextFloat |
|||
| .counts[ (.nextFloat * 5) | floor ] += 1 ) |
|||
| "\nThe counts for \($n) repetitions are:", |
|||
(range(0;5) as $i | " \($i) : \(.counts[$i] // 0)"); |
|||
task1, |
|||
task2(100000) |
|||
</syntaxhighlight> |
|||
{{output}} |
|||
<pre> |
|||
1459213977 |
|||
2827710106 |
|||
4245671317 |
|||
3877608661 |
|||
2595287583 |
|||
The counts for 100000 repetitions are: |
|||
0 : 20002 |
0 : 20002 |
||
1 : 20060 |
1 : 20060 |
||
Line 456: | Line 905: | ||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
< |
<syntaxhighlight lang="julia">const a1 = [0, 1403580, -810728] |
||
const m1 = 2^32 - 209 |
const m1 = 2^32 - 209 |
||
const a2 = [527612, 0, -1370589] |
const a2 = [527612, 0, -1370589] |
||
Line 490: | Line 939: | ||
end |
end |
||
foreach(p -> print(p[1], ": ", p[2], " "), enumerate(hist)) |
foreach(p -> print(p[1], ": ", p[2], " "), enumerate(hist)) |
||
</ |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
1459213977 |
1459213977 |
||
Line 499: | Line 948: | ||
1: 20002 2: 20060 3: 19948 4: 20059 5: 19931 |
1: 20002 2: 20060 3: 19948 4: 20059 5: 19931 |
||
</pre> |
</pre> |
||
=={{header|Kotlin}}== |
|||
{{trans|C++}} |
|||
<syntaxhighlight lang="scala">import kotlin.math.floor |
|||
fun mod(x: Long, y: Long): Long { |
|||
val m = x % y |
|||
return if (m < 0) { |
|||
if (y < 0) { |
|||
m - y |
|||
} else { |
|||
m + y |
|||
} |
|||
} else m |
|||
} |
|||
class RNG { |
|||
// first generator |
|||
private val a1 = arrayOf(0L, 1403580L, -810728L) |
|||
private val m1 = (1L shl 32) - 209 |
|||
private var x1 = arrayOf(0L, 0L, 0L) |
|||
// second generator |
|||
private val a2 = arrayOf(527612L, 0L, -1370589L) |
|||
private val m2 = (1L shl 32) - 22853 |
|||
private var x2 = arrayOf(0L, 0L, 0L) |
|||
private val d = m1 + 1 |
|||
fun seed(state: Long) { |
|||
x1 = arrayOf(state, 0, 0) |
|||
x2 = arrayOf(state, 0, 0) |
|||
} |
|||
fun nextInt(): Long { |
|||
val x1i = mod(a1[0] * x1[0] + a1[1] * x1[1] + a1[2] * x1[2], m1) |
|||
val x2i = mod(a2[0] * x2[0] + a2[1] * x2[1] + a2[2] * x2[2], m2) |
|||
val z = mod(x1i - x2i, m1) |
|||
// keep last three values of the first generator |
|||
x1 = arrayOf(x1i, x1[0], x1[1]) |
|||
// keep last three values of the second generator |
|||
x2 = arrayOf(x2i, x2[0], x2[1]) |
|||
return z + 1 |
|||
} |
|||
fun nextFloat(): Double { |
|||
return nextInt().toDouble() / d |
|||
} |
|||
} |
|||
fun main() { |
|||
val rng = RNG() |
|||
rng.seed(1234567) |
|||
println(rng.nextInt()) |
|||
println(rng.nextInt()) |
|||
println(rng.nextInt()) |
|||
println(rng.nextInt()) |
|||
println(rng.nextInt()) |
|||
println() |
|||
val counts = IntArray(5) |
|||
rng.seed(987654321) |
|||
for (i in 0 until 100_000) { |
|||
val v = floor((rng.nextFloat() * 5.0)).toInt() |
|||
counts[v]++ |
|||
} |
|||
for (iv in counts.withIndex()) { |
|||
println("${iv.index}: ${iv.value}") |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>1459213977 |
|||
2827710106 |
|||
4245671317 |
|||
3877608661 |
|||
2595287583 |
|||
0: 20002 |
|||
1: 20060 |
|||
2: 19948 |
|||
3: 20059 |
|||
4: 19931</pre> |
|||
=={{header|Nim}}== |
|||
<syntaxhighlight lang="nim">import algorithm, math, sequtils, strutils, tables |
|||
const |
|||
# First generator. |
|||
a1 = [int64 0, 1403580, -810728] |
|||
m1: int64 = 2^32 - 209 |
|||
# Second generator. |
|||
a2 = [int64 527612, 0, -1370589] |
|||
m2: int64 = 2^32 - 22853 |
|||
d = m1 + 1 |
|||
type MRG32k3a = object |
|||
x1: array[3, int64] # List of three last values of gen #1. |
|||
x2: array[3, int64] # List of three last values of gen #2. |
|||
func seed(gen: var MRG32k3a; seedState: int64) = |
|||
assert seedState in 1..<d |
|||
gen.x1 = [seedState, 0, 0] |
|||
gen.x2 = [seedState, 0, 0] |
|||
func nextInt(gen: var MRG32k3a): int64 = |
|||
let x1i = floormod(a1[0] * gen.x1[0] + a1[1] * gen.x1[1] + a1[2] * gen.x1[2], m1) |
|||
let x2i = floormod(a2[0] * gen.x2[0] + a2[1] * gen.x2[1] + a2[2] * gen.x2[2], m2) |
|||
# In version 1.4, the following two lines doesn't work. |
|||
# gen.x1 = [x1i, gen.x1[0], gen.x1[1]] # Keep last three. |
|||
# gen.x2 = [x2i, gen.x2[0], gen.x2[1]] # Keep last three. |
|||
gen.x1[2] = gen.x1[1]; gen.x1[1] = gen.x1[0]; gen.x1[0] = x1i |
|||
gen.x2[2] = gen.x2[1]; gen.x2[1] = gen.x2[0]; gen.x2[0] = x2i |
|||
result = floormod(x1i - x2i, m1) + 1 |
|||
func nextFloat(gen: var MRG32k3a): float = |
|||
gen.nextInt().float / d.float |
|||
when isMainModule: |
|||
var gen: MRG32k3a |
|||
gen.seed(1234567) |
|||
for _ in 1..5: |
|||
echo gen.nextInt() |
|||
echo "" |
|||
gen.seed(987654321) |
|||
var counts: CountTable[int] |
|||
for _ in 1..100_000: |
|||
counts.inc int(gen.nextFloat() * 5) |
|||
echo sorted(toSeq(counts.pairs)).mapIt($it[0] & ": " & $it[1]).join(", ")</syntaxhighlight> |
|||
{{out}} |
|||
<pre>1459213977 |
|||
2827710106 |
|||
4245671317 |
|||
3877608661 |
|||
2595287583 |
|||
0: 20002, 1: 20060, 2: 19948, 3: 20059, 4: 19931</pre> |
|||
=={{header|Pari/GP}}== |
|||
Pretty straightforward translation from the directions. Used column/vector multiplication (essentially he dot product) instead of the more tedious form given in the definition of x1i and x2i; rationals (t_FRAC) used in place of floating-point since GP lacks floating-point. |
|||
<syntaxhighlight lang="parigp">a1 = [0, 1403580, -810728]; |
|||
m1 = 2^32-209; |
|||
a2 = [527612, 0, -1370589]; |
|||
m2 = 2^32-22853; |
|||
d = m1+1; |
|||
seed(s)=x1=x2=[s,0,0]; |
|||
next_int()= |
|||
{ |
|||
my(x1i=a1*x1~%m1, x2i=a2*x2~%m2); |
|||
x1 = [x1i, x1[1], x1[2]]; |
|||
x2 = [x2i, x2[1], x2[2]]; |
|||
(x1i-x2i)%m1 + 1; |
|||
} |
|||
next_float()=next_int()/d; |
|||
seed(1234567); |
|||
vector(5,i,next_int()) |
|||
seed(987654321); |
|||
v=vector(5); for(i=1,1e5, v[next_float()*5\1+1]++); v</syntaxhighlight> |
|||
{{out}} |
|||
<pre>%1 = [1459213977, 2827710106, 4245671317, 3877608661, 2595287583] |
|||
%2 = [20002, 20060, 19948, 20059, 19931]</pre> |
|||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
< |
<syntaxhighlight lang="perl">use strict; |
||
use warnings; |
use warnings; |
||
use feature 'say'; |
use feature 'say'; |
||
Line 523: | Line 1,141: | ||
sub next_int { |
sub next_int { |
||
my ( |
my ($self) = @_; |
||
unshift @{$$self{x1}}, ($a1[0] * $$self{x1}[0] + $a1[1] * $$self{x1}[1] + $a1[2] * $$self{x1}[2]) % m1; pop @{$$self{x1}}; |
unshift @{$$self{x1}}, ($a1[0] * $$self{x1}[0] + $a1[1] * $$self{x1}[1] + $a1[2] * $$self{x1}[2]) % m1; pop @{$$self{x1}}; |
||
unshift @{$$self{x2}}, ($a2[0] * $$self{x2}[0] + $a2[1] * $$self{x2}[1] + $a2[2] * $$self{x2}[2]) % m2; pop @{$$self{x2}}; |
unshift @{$$self{x2}}, ($a2[0] * $$self{x2}[0] + $a2[1] * $$self{x2}[1] + $a2[2] * $$self{x2}[2]) % m2; pop @{$$self{x2}}; |
||
Line 529: | Line 1,147: | ||
} |
} |
||
sub next_float { |
sub next_float { $_[0]->next_int / (m1 + 1) } |
||
} |
} |
||
say 'Seed: 1234567, first 5 values:'; |
say 'Seed: 1234567, first 5 values:'; |
||
my $rng = MRG32k3a->new( seed => 1234567 ); |
my $rng = MRG32k3a->new( seed => 1234567 ); |
||
say |
say $rng->next_int for 1..5; |
||
my %h; |
my %h; |
||
say "\nSeed: 987654321, values histogram:"; |
say "\nSeed: 987654321, values histogram:"; |
||
$rng = MRG32k3a->new( seed => 987654321 ); |
$rng = MRG32k3a->new( seed => 987654321 ); |
||
$h{int 5 * |
$h{int 5 * $rng->next_float}++ for 1..100_000; |
||
say "$_ $h{$_}" for sort keys %h;</ |
say "$_ $h{$_}" for sort keys %h;</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Seed: 1234567, first 5 values: |
<pre>Seed: 1234567, first 5 values: |
||
Line 557: | Line 1,175: | ||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<lang Phix>constant |
|||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
|||
-- First generator |
|||
<span style="color: #008080;">constant</span> |
|||
a1 = {0, 1403580, -810728}, |
|||
<span style="color: #000080;font-style:italic;">-- First generator</span> |
|||
m1 = power(2,32) - 209, |
|||
<span style="color: #000000;">a1</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1403580</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">810728</span><span style="color: #0000FF;">},</span> |
|||
-- Second Generator |
|||
<span style="color: #000000;">m1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">32</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">209</span><span style="color: #0000FF;">,</span> |
|||
a2 = {527612, 0, -1370589}, |
|||
<span style="color: #000080;font-style:italic;">-- Second Generator</span> |
|||
m2 = power(2,32) - 22853, |
|||
<span style="color: #000000;">a2</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">527612</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1370589</span><span style="color: #0000FF;">},</span> |
|||
d = m1 + 1 |
|||
<span style="color: #000000;">m2</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">32</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">22853</span><span style="color: #0000FF;">,</span> |
|||
<span style="color: #000000;">d</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">m1</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span> |
|||
sequence x1 = {0, 0, 0}, /* list of three last values of gen #1 */ |
|||
x2 = {0, 0, 0} /* list of three last values of gen #2 */ |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">x1</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},</span> <span style="color: #000080;font-style:italic;">/* list of three last values of gen #1 */</span> |
|||
<span style="color: #000000;">x2</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">}</span> <span style="color: #000080;font-style:italic;">/* list of three last values of gen #2 */</span> |
|||
procedure seed(integer seed_state) |
|||
assert(seed_state>0 and seed_state<d) |
|||
<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;">seed_state</span><span style="color: #0000FF;">)</span> |
|||
x1 = {seed_state, 0, 0} |
|||
<span style="color: #7060A8;">assert</span><span style="color: #0000FF;">(</span><span style="color: #000000;">seed_state</span><span style="color: #0000FF;">></span><span style="color: #000000;">0</span> <span style="color: #008080;">and</span> <span style="color: #000000;">seed_state</span><span style="color: #0000FF;"><</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)</span> |
|||
x2 = {seed_state, 0, 0} |
|||
<span style="color: #000000;">x1</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">seed_state</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">}</span> |
|||
end procedure |
|||
<span style="color: #000000;">x2</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">seed_state</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">}</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> |
|||
function next_int() |
|||
atom x1i = mod(a1[1]*x1[1] + a1[2]*x1[2] + a1[3]*x1[3],m1), |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">next_int</span><span style="color: #0000FF;">()</span> |
|||
x2i = mod(a2[1]*x2[1] + a2[2]*x2[2] + a2[3]*x2[3],m2) |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">x1i</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mod</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a1</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">x1</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">a1</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">x1</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">a1</span><span style="color: #0000FF;">[</span><span style="color: #000000;">3</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">x1</span><span style="color: #0000FF;">[</span><span style="color: #000000;">3</span><span style="color: #0000FF;">],</span><span style="color: #000000;">m1</span><span style="color: #0000FF;">),</span> |
|||
x1 = {x1i, x1[1], x1[2]} /* Keep last three */ |
|||
<span style="color: #000000;">x2i</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mod</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a2</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">x2</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">a2</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">x2</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">a2</span><span style="color: #0000FF;">[</span><span style="color: #000000;">3</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">x2</span><span style="color: #0000FF;">[</span><span style="color: #000000;">3</span><span style="color: #0000FF;">],</span><span style="color: #000000;">m2</span><span style="color: #0000FF;">)</span> |
|||
x2 = {x2i, x2[1], x2[2]} /* Keep last three */ |
|||
<span style="color: #000000;">x1</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x1i</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">x1</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">x1</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]}</span> <span style="color: #000080;font-style:italic;">/* Keep last three */</span> |
|||
atom z = mod(x1i-x2i,m1), |
|||
<span style="color: #000000;">x2</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x2i</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">x2</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">x2</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]}</span> <span style="color: #000080;font-style:italic;">/* Keep last three */</span> |
|||
answer = (z + 1) |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">z</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mod</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x1i</span><span style="color: #0000FF;">-</span><span style="color: #000000;">x2i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m1</span><span style="color: #0000FF;">),</span> |
|||
return answer |
|||
<span style="color: #000000;">answer</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">z</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
end function |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">answer</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
function next_float() |
|||
return next_int() / d |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">next_float</span><span style="color: #0000FF;">()</span> |
|||
end function |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">next_int</span><span style="color: #0000FF;">()</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">d</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
seed(1234567) |
|||
for i=1 to 5 do |
|||
<span style="color: #000000;">seed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1234567</span><span style="color: #0000FF;">)</span> |
|||
printf(1,"%d\n",next_int()) |
|||
<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> |
|||
end for |
|||
<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;">"%d\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">next_int</span><span style="color: #0000FF;">())</span> |
|||
seed(987654321) |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
sequence r = repeat(0,5) |
|||
<span style="color: #000000;">seed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">987654321</span><span style="color: #0000FF;">)</span> |
|||
for i=1 to 100_000 do |
|||
<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> |
|||
r[floor(next_float()*5)+1] += 1 |
|||
<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;">100_000</span> <span style="color: #008080;">do</span> |
|||
end for |
|||
<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</lang> |
|||
<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> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #0000FF;">?</span><span style="color: #000000;">r</span> |
|||
<!--</syntaxhighlight>--> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 607: | Line 1,229: | ||
2595287583 |
2595287583 |
||
{20002,20060,19948,20059,19931} |
{20002,20060,19948,20059,19931} |
||
</pre> |
</pre> |
||
=={{header|Python}}== |
=={{header|Python}}== |
||
< |
<syntaxhighlight lang="python"># Constants |
||
a1 = [0, 1403580, -810728] |
a1 = [0, 1403580, -810728] |
||
m1 = 2**32 - 209 |
m1 = 2**32 - 209 |
||
Line 656: | Line 1,278: | ||
for i in range(100_000): |
for i in range(100_000): |
||
hist[int(random_gen.next_float() *5)] += 1 |
hist[int(random_gen.next_float() *5)] += 1 |
||
print(hist)</ |
print(hist)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 672: | Line 1,294: | ||
All constants are encapsulated within the class. |
All constants are encapsulated within the class. |
||
<lang |
<syntaxhighlight lang="raku" line>class MRG32k3a { |
||
has @!x1; |
has @!x1; |
||
has @!x2; |
has @!x2; |
||
Line 708: | Line 1,330: | ||
say "\nSeed: default; first five Int values:"; |
say "\nSeed: default; first five Int values:"; |
||
$rng = MRG32k3a.new; |
$rng = MRG32k3a.new; |
||
.say for $rng.next-int xx 5;</ |
.say for $rng.next-int xx 5;</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Seed: 1234567; first five Int values: |
<pre>Seed: 1234567; first five Int values: |
||
Line 729: | Line 1,351: | ||
=={{header|Ruby}}== |
=={{header|Ruby}}== |
||
{{trans|C}} |
{{trans|C}} |
||
< |
<syntaxhighlight lang="ruby">def mod(x, y) |
||
m = x % y |
m = x % y |
||
if m < 0 then |
if m < 0 then |
||
Line 796: | Line 1,418: | ||
counts.each_with_index { |v,i| |
counts.each_with_index { |v,i| |
||
print i, ": ", v, "\n" |
print i, ": ", v, "\n" |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>1459213977 |
<pre>1459213977 |
||
Line 809: | Line 1,431: | ||
3: 20059 |
3: 20059 |
||
4: 19931</pre> |
4: 19931</pre> |
||
===Mimicking the Pseudo-code=== |
|||
<syntaxhighlight lang="ruby"># Constants |
|||
# First generator |
|||
A1 = [0, 1403580, -810728] |
|||
M1 = 2**32 - 209 |
|||
# Second Generator |
|||
A2 = [527612, 0, -1370589] |
|||
M2 = 2**32 - 22853 |
|||
D = M1 + 1 |
|||
class MRG32k3a |
|||
def seed(seed_state) |
|||
raise ArgumentError unless seed_state.between?(0, D) |
|||
@x1 = [seed_state, 0, 0] |
|||
@x2 = [seed_state, 0, 0] |
|||
end |
|||
def next_int |
|||
x1i = (A1[0]*@x1[0] + A1[1]*@x1[1] + A1[2]*@x1[2]).modulo M1 |
|||
x2i = (A2[0]*@x2[0] + A2[1]*@x2[1] + A2[2]*@x2[2]).modulo M2 |
|||
@x1 = [x1i, @x1[0], @x1[1]] # Keep last three |
|||
@x2 = [x2i, @x2[0], @x2[1]] # Keep last three |
|||
z = (x1i - x2i) % M1 |
|||
return z + 1 |
|||
end |
|||
def next_float |
|||
next_int.to_f / D |
|||
end |
|||
end |
|||
random_gen = MRG32k3a.new |
|||
random_gen.seed(1234567) |
|||
5.times{ puts random_gen.next_int} |
|||
random_gen = MRG32k3a.new |
|||
random_gen.seed(987654321) |
|||
p 100_000.times.map{(random_gen.next_float() * 5).floor}.tally.sort.to_h |
|||
</syntaxhighlight> |
|||
=={{header|Sidef}}== |
|||
{{trans|Perl}} |
|||
<syntaxhighlight lang="ruby">class MRG32k3a(seed) { |
|||
define( |
|||
m1 = (2**32 - 209) |
|||
m2 = (2**32 - 22853) |
|||
) |
|||
define( |
|||
a1 = %n< 0 1403580 -810728> |
|||
a2 = %n<527612 0 -1370589> |
|||
) |
|||
has x1 = [seed, 0, 0] |
|||
has x2 = x1.clone |
|||
method next_int { |
|||
x1.unshift(a1.map_kv {|k,v| v * x1[k] }.sum % m1); x1.pop |
|||
x2.unshift(a2.map_kv {|k,v| v * x2[k] }.sum % m2); x2.pop |
|||
(x1[0] - x2[0]) % (m1 + 1) |
|||
} |
|||
method next_float { self.next_int / (m1 + 1) -> float } |
|||
} |
|||
say "Seed: 1234567, first 5 values:" |
|||
var rng = MRG32k3a(seed: 1234567) |
|||
5.of { rng.next_int }.each { .say } |
|||
say "\nSeed: 987654321, values histogram:"; |
|||
var rng = MRG32k3a(seed: 987654321) |
|||
var freq = 100_000.of { rng.next_float * 5 -> int }.freq |
|||
freq.sort.each_2d {|k,v| say "#{k} #{v}" }</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Seed: 1234567, first 5 values: |
|||
1459213977 |
|||
2827710106 |
|||
4245671317 |
|||
3877608661 |
|||
2595287583 |
|||
Seed: 987654321, values histogram: |
|||
0 20002 |
|||
1 20060 |
|||
2 19948 |
|||
3 20059 |
|||
4 19931 |
|||
</pre> |
|||
=={{header|uBasic/4tH}}== |
|||
{{works with|v3.64}} |
|||
{{trans|C}} |
|||
Since uBasic/4tH has no floating point support, only the integer part of the task can be implemented. |
|||
<syntaxhighlight lang="text">@(0) = 0 ' First generator |
|||
@(1) = 1403580 |
|||
@(2) = -810728 |
|||
m = SHL(1, 32) - 209 |
|||
@(3) = 527612 ' Second generator |
|||
@(4) = 0 |
|||
@(5) = -1370589 |
|||
n = SHL(1, 32) - 22853 |
|||
d = SHL(1, 32) - 209 + 1 ' m + 1 |
|||
Proc _Seed(1234567) |
|||
Print FUNC(_NextInt) |
|||
Print FUNC(_NextInt) |
|||
Print FUNC(_NextInt) |
|||
Print FUNC(_NextInt) |
|||
Print FUNC(_NextInt) |
|||
Print |
|||
End |
|||
_Mod Param(2) |
|||
Local(1) |
|||
c@ = a@ % b@ |
|||
If c@ < 0 Then |
|||
If b@ < 0 Then |
|||
Return (c@-b@) |
|||
Else |
|||
Return (c@+b@) |
|||
Endif |
|||
EndIf |
|||
Return (c@) |
|||
_Seed Param(1) ' seed the PRNG |
|||
@(6) = a@ |
|||
@(7) = 0 |
|||
@(8) = 0 |
|||
@(9) = a@ |
|||
@(10) = 0 |
|||
@(11) = 0 |
|||
Return |
|||
_NextInt ' get the next random integer value |
|||
Local(3) |
|||
a@ = FUNC(_Mod((@(0) * @(6) + @(1) * @(7) + @(2) * @(8)), m)) |
|||
b@ = FUNC(_Mod((@(3) * @(9) + @(4) * @(10) + @(5) * @(11)), n)) |
|||
c@ = FUNC(_Mod(a@ - b@, m)) |
|||
' keep last three values of the first generator |
|||
@(8) = @(7) |
|||
@(7) = @(6) |
|||
@(6) = a@ |
|||
' keep last three values of the second generator |
|||
@(11) = @(10) |
|||
@(10) = @(9) |
|||
@(9) = b@ |
|||
Return (c@ + 1)</syntaxhighlight> |
|||
{{out}} |
|||
<pre>1459213977 |
|||
2827710106 |
|||
4245671317 |
|||
3877608661 |
|||
2595287583 |
|||
0 OK, 0:398 |
|||
</pre> |
|||
=={{header|Wren}}== |
=={{header|Wren}}== |
||
{{trans|Python}} |
{{trans|Python}} |
||
< |
<syntaxhighlight lang="wren">// constants |
||
var A1 = [0, 1403580, -810728] |
var A1 = [0, 1403580, -810728] |
||
var M1 = 2.pow(32) - 209 |
var M1 = 2.pow(32) - 209 |
||
Line 861: | Line 1,652: | ||
} |
} |
||
System.print("\nThe counts for 100,000 repetitions are:") |
System.print("\nThe counts for 100,000 repetitions are:") |
||
for (i in 0..4) System.print(" %(i) : %(counts[i])")</ |
for (i in 0..4) System.print(" %(i) : %(counts[i])")</syntaxhighlight> |
||
{{out}} |
{{out}} |
Latest revision as of 02:42, 9 June 2024
You are encouraged to solve this task according to the task description, using any language you may know.
- MRG32k3a Combined recursive generator (pseudo-code)
/* Constants */ /* First generator */ a1 = [0, 1403580, -810728] m1 = 2**32 - 209 /* Second Generator */ a2 = [527612, 0, -1370589] m2 = 2**32 - 22853 d = m1 + 1 class MRG32k3a x1 = [0, 0, 0] /* list of three last values of gen #1 */ x2 = [0, 0, 0] /* list of three last values of gen #2 */ method seed(u64 seed_state) assert seed_state in range >0 and < d x1 = [seed_state, 0, 0] x2 = [seed_state, 0, 0] end method method next_int() x1i = (a1[0]*x1[0] + a1[1]*x1[1] + a1[2]*x1[2]) mod m1 x2i = (a2[0]*x2[0] + a2[1]*x2[1] + a2[2]*x2[2]) mod m2 x1 = [x1i, x1[0], x1[1]] /* Keep last three */ x2 = [x2i, x2[0], x2[1]] /* Keep last three */ z = (x1i - x2i) % m1 answer = (z + 1) return answer end method method next_float(): return float next_int() / d end method end class
- MRG32k3a Use:
random_gen = instance MRG32k3a random_gen.seed(1234567) print(random_gen.next_int()) /* 1459213977 */ print(random_gen.next_int()) /* 2827710106 */ print(random_gen.next_int()) /* 4245671317 */ print(random_gen.next_int()) /* 3877608661 */ print(random_gen.next_int()) /* 2595287583 */
- Task
- Generate a class/set of functions that generates pseudo-random
numbers as shown above.
- Show that the first five integers generated 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: 20002, 1: 20060, 2: 19948, 3: 20059, 4: 19931
- Show your output here, on this page.
11l
V a1 = [Int64(0), 1403580, -810728]
V m1 = Int64(2) ^ 32 - 209
V a2 = [Int64(527612), 0, -1370589]
V m2 = Int64(2) ^ 32 - 22853
V d = m1 + 1
T MRG32k3a
[Int64] x1, x2
F (seed_state = 123)
.seed(seed_state)
F seed(Int64 seed_state)
assert(seed_state C Int64(0) <.< :d, ‘Out of Range 0 x < #.’.format(:d))
.x1 = [Int64(seed_state), 0, 0]
.x2 = [Int64(seed_state), 0, 0]
F next_int()
‘return random int in range 0..d’
V x1i = (sum(zip(:a1, .x1).map((aa, xx) -> aa * xx)) % :m1 + :m1) % :m1
V x2i = (sum(zip(:a2, .x2).map((aa, xx) -> aa * xx)) % :m2 + :m2) % :m2
.x1 = [x1i] [+] .x1[0.<2]
.x2 = [x2i] [+] .x2[0.<2]
V z = ((x1i - x2i) % :m1 + :m1) % :m1
R z + 1
F next_float()
‘return random float between 0 and 1’
R Float(.next_int()) / :d
V random_gen = MRG32k3a()
random_gen.seed(1234567)
L 5
print(random_gen.next_int())
random_gen.seed(987654321)
V hist = Dict(0.<5, i -> (i, 0))
L 100'000
hist[Int(random_gen.next_float() * 5)]++
print(hist)
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 [0 = 20002, 1 = 20060, 2 = 19948, 3 = 20059, 4 = 19931]
Ada
package MRG32KA is
type I64 is range -2**63..2**63 - 1;
m1 : constant I64 := 2**32 - 209;
m2 : constant I64 := 2**32 - 22853;
subtype state_value is I64 range 1..m1;
procedure Seed (seed_state : state_value);
function Next_Int return I64;
function Next_Float return Long_Float;
end MRG32KA;
package body MRG32KA is
type Data_Array is array (0..2) of I64;
d : constant I64 := m1 + 1;
----------------
-- Generators --
----------------
a1 : Data_Array := (0, 1403580, -810728);
a2 : Data_Array := (527612, 0, -1370589);
x1 : Data_Array := (0, 0, 0);
x2 : Data_Array := (0, 0, 0);
----------
-- Seed --
----------
procedure Seed (seed_state : state_value) is
begin
x1 := (seed_state, 0, 0);
x2 := (seed_state, 0, 0);
end Seed;
--------------
-- Next_Int --
--------------
function Next_Int return I64 is
x1i : i64;
x2i : I64;
z : I64;
answer : I64;
begin
x1i := (a1(0) * x1(0) + a1(1) * x1(1) + a1(2) * x1(2)) mod m1;
x2i := (a2(0) * x2(0) + a2(1) * x2(1) + a2(2) * x2(2)) mod m2;
x1 := (x1i, x1(0), x1(1));
x2 := (x2i, x2(0), x2(1));
z := (x1i - x2i) mod m1;
answer := z + 1;
return answer;
end Next_Int;
----------------
-- Next_Float --
----------------
function Next_Float return Long_Float is
begin
return Long_float(Next_Int) / Long_Float(d);
end Next_Float;
end MRG32KA;
with Ada.Text_IO; use Ada.Text_IO;
with mrg32ka; use mrg32ka;
procedure Main is
counts : array(0..4) of Natural := (Others => 0);
J : Natural;
begin
seed(1234567);
for I in 1..5 loop
Put_Line(I64'Image(Next_Int));
end loop;
New_Line;
seed(987654321);
for I in 1..100_000 loop
J := Natural(Long_Float'Floor(Next_Float * 5.0));
Counts(J) := Counts(J) + 1;
end loop;
for I in Counts'Range loop
Put(I'Image & " :" & Counts(I)'Image);
end loop;
end Main;
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0 : 20002 1 : 20060 2 : 19948 3 : 20059 4 : 19931
C
#include <math.h>
#include <stdio.h>
#include <stdint.h>
int64_t mod(int64_t x, int64_t y) {
int64_t m = x % y;
if (m < 0) {
if (y < 0) {
return m - y;
} else {
return m + y;
}
}
return m;
}
// Constants
// First generator
const static int64_t a1[3] = { 0, 1403580, -810728 };
const static int64_t m1 = (1LL << 32) - 209;
// Second generator
const static int64_t a2[3] = { 527612, 0, -1370589 };
const static int64_t m2 = (1LL << 32) - 22853;
const static int64_t d = (1LL << 32) - 209 + 1; // m1 + 1
// the last three values of the first generator
static int64_t x1[3];
// the last three values of the second generator
static int64_t x2[3];
void seed(int64_t seed_state) {
x1[0] = seed_state;
x1[1] = 0;
x1[2] = 0;
x2[0] = seed_state;
x2[1] = 0;
x2[2] = 0;
}
int64_t next_int() {
int64_t x1i = mod((a1[0] * x1[0] + a1[1] * x1[1] + a1[2] * x1[2]), m1);
int64_t x2i = mod((a2[0] * x2[0] + a2[1] * x2[1] + a2[2] * x2[2]), m2);
int64_t z = mod(x1i - x2i, m1);
// keep last three values of the first generator
x1[2] = x1[1];
x1[1] = x1[0];
x1[0] = x1i;
// keep last three values of the second generator
x2[2] = x2[1];
x2[1] = x2[0];
x2[0] = x2i;
return z + 1;
}
double next_float() {
return (double)next_int() / d;
}
int main() {
int counts[5] = { 0, 0, 0, 0, 0 };
int i;
seed(1234567);
printf("%lld\n", next_int());
printf("%lld\n", next_int());
printf("%lld\n", next_int());
printf("%lld\n", next_int());
printf("%lld\n", next_int());
printf("\n");
seed(987654321);
for (i = 0; i < 100000; i++) {
int64_t value = floor(next_float() * 5);
counts[value]++;
}
for (i = 0; i < 5; i++) {
printf("%d: %d\n", i, counts[i]);
}
return 0;
}
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0: 20002 1: 20060 2: 19948 3: 20059 4: 19931
C++
#include <array>
#include <iostream>
int64_t mod(int64_t x, int64_t y) {
int64_t m = x % y;
if (m < 0) {
if (y < 0) {
return m - y;
} else {
return m + y;
}
}
return m;
}
class RNG {
private:
// First generator
const std::array<int64_t, 3> a1{ 0, 1403580, -810728 };
const int64_t m1 = (1LL << 32) - 209;
std::array<int64_t, 3> x1;
// Second generator
const std::array<int64_t, 3> a2{ 527612, 0, -1370589 };
const int64_t m2 = (1LL << 32) - 22853;
std::array<int64_t, 3> x2;
// other
const int64_t d = (1LL << 32) - 209 + 1; // m1 + 1
public:
void seed(int64_t state) {
x1 = { state, 0, 0 };
x2 = { state, 0, 0 };
}
int64_t next_int() {
int64_t x1i = mod((a1[0] * x1[0] + a1[1] * x1[1] + a1[2] * x1[2]), m1);
int64_t x2i = mod((a2[0] * x2[0] + a2[1] * x2[1] + a2[2] * x2[2]), m2);
int64_t z = mod(x1i - x2i, m1);
// keep last three values of the first generator
x1 = { x1i, x1[0], x1[1] };
// keep last three values of the second generator
x2 = { x2i, x2[0], x2[1] };
return z + 1;
}
double next_float() {
return static_cast<double>(next_int()) / d;
}
};
int main() {
RNG rng;
rng.seed(1234567);
std::cout << rng.next_int() << '\n';
std::cout << rng.next_int() << '\n';
std::cout << rng.next_int() << '\n';
std::cout << rng.next_int() << '\n';
std::cout << rng.next_int() << '\n';
std::cout << '\n';
std::array<int, 5> counts{ 0, 0, 0, 0, 0 };
rng.seed(987654321);
for (size_t i = 0; i < 100000; i++) {
auto value = floor(rng.next_float() * 5.0);
counts[value]++;
}
for (size_t i = 0; i < counts.size(); i++) {
std::cout << i << ": " << counts[i] << '\n';
}
return 0;
}
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0: 20002 1: 20060 2: 19948 3: 20059 4: 19931
D
import std.math;
import std.stdio;
long mod(long x, long y) {
long m = x % y;
if (m < 0) {
if (y < 0) {
return m - y;
} else {
return m + y;
}
}
return m;
}
class RNG {
private:
// First generator
immutable(long []) a1 = [0, 1403580, -810728];
immutable long m1 = (1L << 32) - 209;
long[3] x1;
// Second generator
immutable(long []) a2 = [527612, 0, -1370589];
immutable long m2 = (1L << 32) - 22853;
long[3] x2;
// other
immutable long d = m1 + 1;
public:
void seed(long state) {
x1 = [state, 0, 0];
x2 = [state, 0, 0];
}
long next_int() {
long x1i = mod((a1[0] * x1[0] + a1[1] * x1[1] + a1[2] * x1[2]), m1);
long x2i = mod((a2[0] * x2[0] + a2[1] * x2[1] + a2[2] * x2[2]), m2);
long z = mod(x1i - x2i, m1);
// keep the last three values of the first generator
x1 = [x1i, x1[0], x1[1]];
// keep the last three values of the second generator
x2 = [x2i, x2[0], x2[1]];
return z + 1;
}
double next_float() {
return cast(double) next_int() / d;
}
}
void main() {
auto rng = new RNG();
rng.seed(1234567);
writeln(rng.next_int);
writeln(rng.next_int);
writeln(rng.next_int);
writeln(rng.next_int);
writeln(rng.next_int);
writeln;
int[5] counts;
rng.seed(987654321);
foreach (i; 0 .. 100_000) {
auto value = cast(int) floor(rng.next_float * 5.0);
counts[value]++;
}
foreach (i,v; counts) {
writeln(i, ": ", v);
}
}
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0: 20002 1: 20060 2: 19948 3: 20059 4: 19931
Factor
USING: arrays kernel math math.order math.statistics
math.vectors prettyprint sequences ;
CONSTANT: m1 4294967087
CONSTANT: m2 4294944443
: seed ( n -- seq1 seq2 )
dup 1 m1 between? t assert= 0 0 3array dup ;
: new-state ( seq1 seq2 n -- new-seq )
[ dup ] [ vdot ] [ rem prefix but-last ] tri* ;
: next-state ( a b -- a' b' )
[ { 0 1403580 -810728 } m1 new-state ]
[ { 527612 0 -1370589 } m2 new-state ] bi* ;
: next-int ( a b -- a' b' n )
next-state 2dup [ first ] bi@ - m1 rem 1 + ;
: next-float ( a b -- a' b' x ) next-int m1 1 + /f ;
! Task
1234567 seed 5 [ next-int . ] times 2drop
987654321 seed 100,000 [ next-float 5 * >integer ] replicate
2nip histogram .
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 H{ { 0 20002 } { 1 20060 } { 2 19948 } { 3 20059 } { 4 19931 } }
Forth
6 array (seed) \ holds the seed
6 array (gens) \ holds the generators
\ set up constants
0 (gens) 0 th ! \ 1st generator
1403580 (gens) 1 th !
-810728 (gens) 2 th !
527612 (gens) 3 th ! \ 2nd generator
0 (gens) 4 th !
-1370589 (gens) 5 th !
1 32 lshift 209 - value (m) \ 1st generator constant
1 32 lshift 22853 - value (n) \ 2nd generator constant
( n1 n2 -- n3)
: (mod) tuck mod tuck 0< if abs + ;then drop ;
: (generate) do (seed) i th @ (gens) i th @ * + loop swap (mod) ;
: (reseed) ?do (seed) i th ! loop ; ( n1 n2 n3 limit index --)
: randomize 6 0 do dup i 3 mod if >zero then (seed) i th ! loop drop ;
( n --)
: random ( -- n)
(m) 0 3 0 (generate) (n) 0 6 3 (generate) over over
(seed) 4 th @ (seed) 3 th @ rot 6 3 (reseed)
(seed) 1 th @ (seed) 0 th @ rot 3 0 (reseed) - (m) (mod) 1+
;
include lib/fp1.4th \ simple floating point support
include lib/zenfloor.4th \ for FLOOR
5 array (count) \ setup an array of 5 elements
: test
1234567 randomize
random . cr random . cr random . cr
random . cr random . cr cr \ perform the first test
987654321 randomize (m) 1+ s>f \ set up denominator
100000 0 ?do \ do this 100,000 times
random s>f fover f/ 5 s>f f* floor f>s cells (count) + 1 swap +!
loop fdrop
\ show the results
5 0 ?do i . ." : " (count) i th ? cr loop
;
test
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0 : 20002 1 : 20060 2 : 19948 3 : 20059 4 : 19931
Go
package main
import (
"fmt"
"log"
"math"
)
var a1 = []int64{0, 1403580, -810728}
var a2 = []int64{527612, 0, -1370589}
const m1 = int64((1 << 32) - 209)
const m2 = int64((1 << 32) - 22853)
const d = m1 + 1
// Python style modulus
func mod(x, y int64) int64 {
m := x % y
if m < 0 {
if y < 0 {
return m - y
} else {
return m + y
}
}
return m
}
type MRG32k3a struct{ x1, x2 [3]int64 }
func MRG32k3aNew() *MRG32k3a { return &MRG32k3a{} }
func (mrg *MRG32k3a) seed(seedState int64) {
if seedState <= 0 || seedState >= d {
log.Fatalf("Argument must be in the range [0, %d].\n", d)
}
mrg.x1 = [3]int64{seedState, 0, 0}
mrg.x2 = [3]int64{seedState, 0, 0}
}
func (mrg *MRG32k3a) nextInt() int64 {
x1i := mod(a1[0]*mrg.x1[0]+a1[1]*mrg.x1[1]+a1[2]*mrg.x1[2], m1)
x2i := mod(a2[0]*mrg.x2[0]+a2[1]*mrg.x2[1]+a2[2]*mrg.x2[2], m2)
mrg.x1 = [3]int64{x1i, mrg.x1[0], mrg.x1[1]} /* keep last three */
mrg.x2 = [3]int64{x2i, mrg.x2[0], mrg.x2[1]} /* keep last three */
return mod(x1i-x2i, m1) + 1
}
func (mrg *MRG32k3a) nextFloat() float64 { return float64(mrg.nextInt()) / float64(d) }
func main() {
randomGen := MRG32k3aNew()
randomGen.seed(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:
1459213977 2827710106 4245671317 3877608661 2595287583 The counts for 100,000 repetitions are: 0 : 20002 1 : 20060 2 : 19948 3 : 20059 4 : 19931
Haskell
import Data.List
randoms :: Int -> [Int]
randoms seed = unfoldr go ([seed,0,0],[seed,0,0])
where
go (x1,x2) =
let x1i = sum (zipWith (*) x1 a1) `mod` m1
x2i = sum (zipWith (*) x2 a2) `mod` m2
in Just $ ((x1i - x2i) `mod` m1, (x1i:init x1, x2i:init x2))
a1 = [0, 1403580, -810728]
m1 = 2^32 - 209
a2 = [527612, 0, -1370589]
m2 = 2^32 - 22853
randomsFloat = map ((/ (2^32 - 208)) . fromIntegral) . randoms
*Main> take 5 $ randoms 1234567 [1459213976,2827710105,4245671316,3877608660,2595287582] *Main> let hist = map length . group . sort *Main> hist . take 100000 $ (floor . (*5)) <$> randomsFloat 987654321 [20002,20060,19948,20059,19931]
As a RandomGen instanse
import System.Random
newtype MRG32k3a = MRG32k3a ([Int],[Int])
mkMRG32k3a s = MRG32k3a ([s,0,0],[s,0,0])
instance RandomGen MRG32k3a where
next (MRG32k3a (x1,x2)) =
let x1i = sum (zipWith (*) x1 a1) `mod` m1
x2i = sum (zipWith (*) x2 a2) `mod` m2
in ((x1i - x2i) `mod` m1, MRG32k3a (x1i:init x1, x2i:init x2))
where
a1 = [0, 1403580, -810728]
m1 = 2^32 - 209
a2 = [527612, 0, -1370589]
m2 = 2^32 - 22853
split _ = error "MRG32k3a is not splittable"
In this case the sequence or numbers differs from direct unfolding, due to internal uniform shuffling.
*Main> take 5 $ randoms (mkMRG32k3a 1234567) [2827710105,3877608660,3642754129,1259674122,3002249941] *Main> let hist = map length . group . sort *Main> hist . take 100000 $ (floor . (*5)) <$> (randoms (mkMRG32k3a 987654321) :: [Float]) [20015,19789,20024,20133,20039]
Java
public class App {
private static long mod(long x, long y) {
long m = x % y;
if (m < 0) {
if (y < 0) {
return m - y;
} else {
return m + y;
}
}
return m;
}
public static class RNG {
// first generator
private final long[] a1 = {0, 1403580, -810728};
private static final long m1 = (1L << 32) - 209;
private long[] x1;
// second generator
private final long[] a2 = {527612, 0, -1370589};
private static final long m2 = (1L << 32) - 22853;
private long[] x2;
// other
private static final long d = m1 + 1;
public void seed(long state) {
x1 = new long[]{state, 0, 0};
x2 = new long[]{state, 0, 0};
}
public long nextInt() {
long x1i = mod(a1[0] * x1[0] + a1[1] * x1[1] + a1[2] * x1[2], m1);
long x2i = mod(a2[0] * x2[0] + a2[1] * x2[1] + a2[2] * x2[2], m2);
long z = mod(x1i - x2i, m1);
// keep the last three values of the first generator
x1 = new long[]{x1i, x1[0], x1[1]};
// keep the last three values of the second generator
x2 = new long[]{x2i, x2[0], x2[1]};
return z + 1;
}
public double nextFloat() {
return 1.0 * nextInt() / d;
}
}
public static void main(String[] args) {
RNG rng = new RNG();
rng.seed(1234567);
System.out.println(rng.nextInt());
System.out.println(rng.nextInt());
System.out.println(rng.nextInt());
System.out.println(rng.nextInt());
System.out.println(rng.nextInt());
System.out.println();
int[] counts = {0, 0, 0, 0, 0};
rng.seed(987654321);
for (int i = 0; i < 100_000; i++) {
int value = (int) Math.floor(rng.nextFloat() * 5.0);
counts[value]++;
}
for (int i = 0; i < counts.length; i++) {
System.out.printf("%d: %d%n", i, counts[i]);
}
}
}
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0: 20002 1: 20060 2: 19948 3: 20059 4: 19931
jq
Adapted from Wren
Works with jq, the C implementation of jq
Works with gojq, the Go implementation of jq
Works with jaq, the Rust implementation of jq
# constants
def A1: [0, 1403580, -810728];
def M1: 4294967087; # pow(2;32) - 209
def A2: [527612, 0, -1370589];
def M2: 4294944443; # pow(2;32) - 22853
def D: M1 + 1;
# Python-style modulus
def Mod($x; $y):
def abs: if . < 0 then - . else . end;
($x % $y) as $m
| if $m < 0 then $m + ($y|abs) else $m end;
def MRG32k3a: {x1: [0, 0, 0], x2: [0, 0, 0] };
def seed($seedState):
if $seedState <= 0 or $seedState >= D
then "Argument of seed must be in the range (0, \(D))" | error
else {x1: [$seedState, 0, 0]} | .x2 = .x1
end ;
# Input: {x1, x2} as for MRG32k31
# Output: {x1, x2, nextInt}
def nextInt:
Mod(A1[0]*.x1[0] + A1[1]*.x1[1] + A1[2]*.x1[2]; M1) as $x1i
| Mod(A2[0]*.x2[0] + A2[1]*.x2[1] + A2[2]*.x2[2]; M2) as $x2i
| .x1 = [$x1i, .x1[0], .x1[1]] # keep last three
| .x2 = [$x2i, .x2[0], .x2[1]] # keep last three
| .nextInt = Mod($x1i - $x2i; M1) + 1 ;
def nextFloat: nextInt | .nextFloat = (.nextInt / D);
def task1:
foreach range(0; 5) as $i (seed(1234567); nextInt ) | .nextInt;
def task2($n):
seed(987654321)
| reduce range(0; $n) as $i (.counts = [range(0;5)|0];
nextFloat
| .counts[ (.nextFloat * 5) | floor ] += 1 )
| "\nThe counts for \($n) repetitions are:",
(range(0;5) as $i | " \($i) : \(.counts[$i] // 0)");
task1,
task2(100000)
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 The counts for 100000 repetitions are: 0 : 20002 1 : 20060 2 : 19948 3 : 20059 4 : 19931
Julia
const a1 = [0, 1403580, -810728]
const m1 = 2^32 - 209
const a2 = [527612, 0, -1370589]
const m2 = 2^32 - 22853
const d = m1 + 1
mutable struct MRG32k3a
x1::Tuple{Int64, Int64, Int64}
x2::Tuple{Int64, Int64, Int64}
MRG32k3a() = new((0, 0, 0), (0, 0, 0))
MRG32k3a(seed_state) = new((seed_state, 0, 0), (seed_state, 0, 0))
end
seed(sd) = begin @assert(0 < sd < d); MRG32k3a(sd) end
function next_int(x::MRG32k3a)
x1i = mod1(a1[1] * x.x1[1] + a1[2] * x.x1[2] + a1[3] * x.x1[3], m1)
x2i = mod1(a2[1] * x.x2[1] + a2[2] * x.x2[2] + a2[3] * x.x2[3], m2)
x.x1 = (x1i, x.x1[1], x.x1[2])
x.x2 = (x2i, x.x2[1], x.x2[2])
return mod1(x1i - x2i, m1) + 1
end
next_float(x::MRG32k3a) = next_int(x) / d
const g1 = seed(1234567)
for _ in 1:5
println(next_int(g1))
end
const g2 = seed(987654321)
hist = fill(0, 5)
for _ in 1:100_000
hist[Int(floor(next_float(g2) * 5)) + 1] += 1
end
foreach(p -> print(p[1], ": ", p[2], " "), enumerate(hist))
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 1: 20002 2: 20060 3: 19948 4: 20059 5: 19931
Kotlin
import kotlin.math.floor
fun mod(x: Long, y: Long): Long {
val m = x % y
return if (m < 0) {
if (y < 0) {
m - y
} else {
m + y
}
} else m
}
class RNG {
// first generator
private val a1 = arrayOf(0L, 1403580L, -810728L)
private val m1 = (1L shl 32) - 209
private var x1 = arrayOf(0L, 0L, 0L)
// second generator
private val a2 = arrayOf(527612L, 0L, -1370589L)
private val m2 = (1L shl 32) - 22853
private var x2 = arrayOf(0L, 0L, 0L)
private val d = m1 + 1
fun seed(state: Long) {
x1 = arrayOf(state, 0, 0)
x2 = arrayOf(state, 0, 0)
}
fun nextInt(): Long {
val x1i = mod(a1[0] * x1[0] + a1[1] * x1[1] + a1[2] * x1[2], m1)
val x2i = mod(a2[0] * x2[0] + a2[1] * x2[1] + a2[2] * x2[2], m2)
val z = mod(x1i - x2i, m1)
// keep last three values of the first generator
x1 = arrayOf(x1i, x1[0], x1[1])
// keep last three values of the second generator
x2 = arrayOf(x2i, x2[0], x2[1])
return z + 1
}
fun nextFloat(): Double {
return nextInt().toDouble() / d
}
}
fun main() {
val rng = RNG()
rng.seed(1234567)
println(rng.nextInt())
println(rng.nextInt())
println(rng.nextInt())
println(rng.nextInt())
println(rng.nextInt())
println()
val counts = IntArray(5)
rng.seed(987654321)
for (i in 0 until 100_000) {
val v = floor((rng.nextFloat() * 5.0)).toInt()
counts[v]++
}
for (iv in counts.withIndex()) {
println("${iv.index}: ${iv.value}")
}
}
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0: 20002 1: 20060 2: 19948 3: 20059 4: 19931
Nim
import algorithm, math, sequtils, strutils, tables
const
# First generator.
a1 = [int64 0, 1403580, -810728]
m1: int64 = 2^32 - 209
# Second generator.
a2 = [int64 527612, 0, -1370589]
m2: int64 = 2^32 - 22853
d = m1 + 1
type MRG32k3a = object
x1: array[3, int64] # List of three last values of gen #1.
x2: array[3, int64] # List of three last values of gen #2.
func seed(gen: var MRG32k3a; seedState: int64) =
assert seedState in 1..<d
gen.x1 = [seedState, 0, 0]
gen.x2 = [seedState, 0, 0]
func nextInt(gen: var MRG32k3a): int64 =
let x1i = floormod(a1[0] * gen.x1[0] + a1[1] * gen.x1[1] + a1[2] * gen.x1[2], m1)
let x2i = floormod(a2[0] * gen.x2[0] + a2[1] * gen.x2[1] + a2[2] * gen.x2[2], m2)
# In version 1.4, the following two lines doesn't work.
# gen.x1 = [x1i, gen.x1[0], gen.x1[1]] # Keep last three.
# gen.x2 = [x2i, gen.x2[0], gen.x2[1]] # Keep last three.
gen.x1[2] = gen.x1[1]; gen.x1[1] = gen.x1[0]; gen.x1[0] = x1i
gen.x2[2] = gen.x2[1]; gen.x2[1] = gen.x2[0]; gen.x2[0] = x2i
result = floormod(x1i - x2i, m1) + 1
func nextFloat(gen: var MRG32k3a): float =
gen.nextInt().float / d.float
when isMainModule:
var gen: MRG32k3a
gen.seed(1234567)
for _ in 1..5:
echo gen.nextInt()
echo ""
gen.seed(987654321)
var counts: CountTable[int]
for _ in 1..100_000:
counts.inc int(gen.nextFloat() * 5)
echo sorted(toSeq(counts.pairs)).mapIt($it[0] & ": " & $it[1]).join(", ")
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0: 20002, 1: 20060, 2: 19948, 3: 20059, 4: 19931
Pari/GP
Pretty straightforward translation from the directions. Used column/vector multiplication (essentially he dot product) instead of the more tedious form given in the definition of x1i and x2i; rationals (t_FRAC) used in place of floating-point since GP lacks floating-point.
a1 = [0, 1403580, -810728];
m1 = 2^32-209;
a2 = [527612, 0, -1370589];
m2 = 2^32-22853;
d = m1+1;
seed(s)=x1=x2=[s,0,0];
next_int()=
{
my(x1i=a1*x1~%m1, x2i=a2*x2~%m2);
x1 = [x1i, x1[1], x1[2]];
x2 = [x2i, x2[1], x2[2]];
(x1i-x2i)%m1 + 1;
}
next_float()=next_int()/d;
seed(1234567);
vector(5,i,next_int())
seed(987654321);
v=vector(5); for(i=1,1e5, v[next_float()*5\1+1]++); v
- Output:
%1 = [1459213977, 2827710106, 4245671317, 3877608661, 2595287583] %2 = [20002, 20060, 19948, 20059, 19931]
Perl
use strict;
use warnings;
use feature 'say';
package MRG32k3a {
use constant {
m1 => 2**32 - 209,
m2 => 2**32 - 22853
};
use Const::Fast;
const my @a1 => < 0 1403580 -810728>;
const my @a2 => <527612 0 -1370589>;
sub new {
my ($class,undef,$seed) = @_;
my @x1 = my @x2 = ($seed, 0, 0);
bless {x1 => \@x1, x2 => \@x2}, $class;
}
sub next_int {
my ($self) = @_;
unshift @{$$self{x1}}, ($a1[0] * $$self{x1}[0] + $a1[1] * $$self{x1}[1] + $a1[2] * $$self{x1}[2]) % m1; pop @{$$self{x1}};
unshift @{$$self{x2}}, ($a2[0] * $$self{x2}[0] + $a2[1] * $$self{x2}[1] + $a2[2] * $$self{x2}[2]) % m2; pop @{$$self{x2}};
($$self{x1}[0] - $$self{x2}[0]) % (m1 + 1)
}
sub next_float { $_[0]->next_int / (m1 + 1) }
}
say 'Seed: 1234567, first 5 values:';
my $rng = MRG32k3a->new( seed => 1234567 );
say $rng->next_int for 1..5;
my %h;
say "\nSeed: 987654321, values histogram:";
$rng = MRG32k3a->new( seed => 987654321 );
$h{int 5 * $rng->next_float}++ for 1..100_000;
say "$_ $h{$_}" for sort keys %h;
- Output:
Seed: 1234567, first 5 values: 1459213977 2827710106 4245671317 3877608661 2595287583 Seed: 987654321, values histogram: 0 20002 1 20060 2 19948 3 20059 4 19931
Phix
with javascript_semantics constant -- First generator a1 = {0, 1403580, -810728}, m1 = power(2,32) - 209, -- Second Generator a2 = {527612, 0, -1370589}, m2 = power(2,32) - 22853, d = m1 + 1 sequence x1 = {0, 0, 0}, /* list of three last values of gen #1 */ x2 = {0, 0, 0} /* list of three last values of gen #2 */ procedure seed(integer seed_state) assert(seed_state>0 and seed_state<d) x1 = {seed_state, 0, 0} x2 = {seed_state, 0, 0} end procedure function next_int() atom x1i = mod(a1[1]*x1[1] + a1[2]*x1[2] + a1[3]*x1[3],m1), x2i = mod(a2[1]*x2[1] + a2[2]*x2[2] + a2[3]*x2[3],m2) x1 = {x1i, x1[1], x1[2]} /* Keep last three */ x2 = {x2i, x2[1], x2[2]} /* Keep last three */ atom z = mod(x1i-x2i,m1), answer = (z + 1) return answer end function function next_float() return next_int() / d 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 100_000 do integer rdx = floor(next_float()*5)+1 r[rdx] += 1 end for ?r
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 {20002,20060,19948,20059,19931}
Python
# Constants
a1 = [0, 1403580, -810728]
m1 = 2**32 - 209
#
a2 = [527612, 0, -1370589]
m2 = 2**32 - 22853
#
d = m1 + 1
class MRG32k3a():
def __init__(self, seed_state=123):
self.seed(seed_state)
def seed(self, seed_state):
assert 0 <seed_state < d, f"Out of Range 0 x < {d}"
self.x1 = [seed_state, 0, 0]
self.x2 = [seed_state, 0, 0]
def next_int(self):
"return random int in range 0..d"
x1i = sum(aa * xx for aa, xx in zip(a1, self.x1)) % m1
x2i = sum(aa * xx for aa, xx in zip(a2, self.x2)) % m2
self.x1 = [x1i] + self.x1[:2]
self.x2 = [x2i] + self.x2[:2]
z = (x1i - x2i) % m1
answer = (z + 1)
return answer
def next_float(self):
"return random float between 0 and 1"
return self.next_int() / d
if __name__ == '__main__':
random_gen = MRG32k3a()
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:
1459213977 2827710106 4245671317 3877608661 2595287583 {0: 20002, 1: 20060, 2: 19948, 3: 20059, 4: 19931}
Raku
All constants are encapsulated within the class.
class MRG32k3a {
has @!x1;
has @!x2;
constant a1 = 0, 1403580, -810728;
constant a2 = 527612, 0, -1370589;
constant m1 = 2**32 - 209;
constant m2 = 2**32 - 22853;
submethod BUILD ( Int :$seed where 0 < * <= m1 = 1 ) { @!x1 = @!x2 = $seed, 0, 0 }
method next-int {
@!x1.unshift: (a1[0] * @!x1[0] + a1[1] * @!x1[1] + a1[2] * @!x1[2]) % m1; @!x1.pop;
@!x2.unshift: (a2[0] * @!x2[0] + a2[1] * @!x2[1] + a2[2] * @!x2[2]) % m2; @!x2.pop;
(@!x1[0] - @!x2[0]) % (m1 + 1)
}
method next-rat { self.next-int / (m1 + 1) }
}
# Test next-int with custom seed
say 'Seed: 1234567; first five Int values:';
my $rng = MRG32k3a.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 = MRG32k3a.new :seed(987654321);
say ( ($rng.next-rat * 5).floor xx 100_000 ).Bag;
# Test next-int with default seed
say "\nSeed: default; first five Int values:";
$rng = MRG32k3a.new;
.say for $rng.next-int xx 5;
- Output:
Seed: 1234567; first five Int values: 1459213977 2827710106 4245671317 3877608661 2595287583 Seed: 987654321; first 1e5 Rat values histogram: Bag(0(20002) 1(20060) 2(19948) 3(20059) 4(19931)) Seed: default; first five Int values: 4294439476 798392476 1012402088 1268414424 3353586348
Ruby
def mod(x, y)
m = x % y
if m < 0 then
if y < 0 then
return m - y
else
return m + y
end
end
return m
end
# Constants
# First generator
A1 = [0, 1403580, -810728]
A1.freeze
M1 = (1 << 32) - 209
# Second generator
A2 = [527612, 0, -1370589]
A2.freeze
M2 = (1 << 32) - 22853
D = M1 + 1
# the last three values of the first generator
$x1 = [0, 0, 0]
# the last three values of the second generator
$x2 = [0, 0, 0]
def seed(seed_state)
$x1 = [seed_state, 0, 0]
$x2 = [seed_state, 0, 0]
end
def next_int()
x1i = mod((A1[0] * $x1[0] + A1[1] * $x1[1] + A1[2] * $x1[2]), M1)
x2i = mod((A2[0] * $x2[0] + A2[1] * $x2[1] + A2[2] * $x2[2]), M2)
z = mod(x1i - x2i, M1)
$x1 = [x1i, $x1[0], $x1[1]]
$x2 = [x2i, $x2[0], $x2[1]]
return z + 1
end
def next_float()
return 1.0 * next_int() / D
end
########################################
seed(1234567)
print next_int(), "\n"
print next_int(), "\n"
print next_int(), "\n"
print next_int(), "\n"
print next_int(), "\n"
print "\n"
counts = [0, 0, 0, 0, 0]
seed(987654321)
for i in 1 .. 100000
value = (next_float() * 5.0).floor
counts[value] = counts[value] + 1
end
counts.each_with_index { |v,i|
print i, ": ", v, "\n"
}
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0: 20002 1: 20060 2: 19948 3: 20059 4: 19931
Mimicking the Pseudo-code
# Constants
# First generator
A1 = [0, 1403580, -810728]
M1 = 2**32 - 209
# Second Generator
A2 = [527612, 0, -1370589]
M2 = 2**32 - 22853
D = M1 + 1
class MRG32k3a
def seed(seed_state)
raise ArgumentError unless seed_state.between?(0, D)
@x1 = [seed_state, 0, 0]
@x2 = [seed_state, 0, 0]
end
def next_int
x1i = (A1[0]*@x1[0] + A1[1]*@x1[1] + A1[2]*@x1[2]).modulo M1
x2i = (A2[0]*@x2[0] + A2[1]*@x2[1] + A2[2]*@x2[2]).modulo M2
@x1 = [x1i, @x1[0], @x1[1]] # Keep last three
@x2 = [x2i, @x2[0], @x2[1]] # Keep last three
z = (x1i - x2i) % M1
return z + 1
end
def next_float
next_int.to_f / D
end
end
random_gen = MRG32k3a.new
random_gen.seed(1234567)
5.times{ puts random_gen.next_int}
random_gen = MRG32k3a.new
random_gen.seed(987654321)
p 100_000.times.map{(random_gen.next_float() * 5).floor}.tally.sort.to_h
Sidef
class MRG32k3a(seed) {
define(
m1 = (2**32 - 209)
m2 = (2**32 - 22853)
)
define(
a1 = %n< 0 1403580 -810728>
a2 = %n<527612 0 -1370589>
)
has x1 = [seed, 0, 0]
has x2 = x1.clone
method next_int {
x1.unshift(a1.map_kv {|k,v| v * x1[k] }.sum % m1); x1.pop
x2.unshift(a2.map_kv {|k,v| v * x2[k] }.sum % m2); x2.pop
(x1[0] - x2[0]) % (m1 + 1)
}
method next_float { self.next_int / (m1 + 1) -> float }
}
say "Seed: 1234567, first 5 values:"
var rng = MRG32k3a(seed: 1234567)
5.of { rng.next_int }.each { .say }
say "\nSeed: 987654321, values histogram:";
var rng = MRG32k3a(seed: 987654321)
var freq = 100_000.of { rng.next_float * 5 -> int }.freq
freq.sort.each_2d {|k,v| say "#{k} #{v}" }
- Output:
Seed: 1234567, first 5 values: 1459213977 2827710106 4245671317 3877608661 2595287583 Seed: 987654321, values histogram: 0 20002 1 20060 2 19948 3 20059 4 19931
uBasic/4tH
Since uBasic/4tH has no floating point support, only the integer part of the task can be implemented.
@(0) = 0 ' First generator
@(1) = 1403580
@(2) = -810728
m = SHL(1, 32) - 209
@(3) = 527612 ' Second generator
@(4) = 0
@(5) = -1370589
n = SHL(1, 32) - 22853
d = SHL(1, 32) - 209 + 1 ' m + 1
Proc _Seed(1234567)
Print FUNC(_NextInt)
Print FUNC(_NextInt)
Print FUNC(_NextInt)
Print FUNC(_NextInt)
Print FUNC(_NextInt)
Print
End
_Mod Param(2)
Local(1)
c@ = a@ % b@
If c@ < 0 Then
If b@ < 0 Then
Return (c@-b@)
Else
Return (c@+b@)
Endif
EndIf
Return (c@)
_Seed Param(1) ' seed the PRNG
@(6) = a@
@(7) = 0
@(8) = 0
@(9) = a@
@(10) = 0
@(11) = 0
Return
_NextInt ' get the next random integer value
Local(3)
a@ = FUNC(_Mod((@(0) * @(6) + @(1) * @(7) + @(2) * @(8)), m))
b@ = FUNC(_Mod((@(3) * @(9) + @(4) * @(10) + @(5) * @(11)), n))
c@ = FUNC(_Mod(a@ - b@, m))
' keep last three values of the first generator
@(8) = @(7)
@(7) = @(6)
@(6) = a@
' keep last three values of the second generator
@(11) = @(10)
@(10) = @(9)
@(9) = b@
Return (c@ + 1)
- Output:
1459213977 2827710106 4245671317 3877608661 2595287583 0 OK, 0:398
Wren
// constants
var A1 = [0, 1403580, -810728]
var M1 = 2.pow(32) - 209
var A2 = [527612, 0, -1370589]
var M2 = 2.pow(32) - 22853
var D = M1 + 1
// Python style modulus
var Mod = Fn.new { |x, y|
var m = x % y
return (m < 0) ? m + y.abs : m
}
class MRG32k3a {
construct new() {
_x1 = [0, 0, 0]
_x2 = [0, 0, 0]
}
seed(seedState) {
if (seedState <= 0 || seedState >= D) {
Fiber.abort("Argument must be in the range [0, %(D)].")
}
_x1 = [seedState, 0, 0]
_x2 = [seedState, 0, 0]
}
nextInt {
var x1i = Mod.call(A1[0]*_x1[0] + A1[1]*_x1[1] + A1[2]*_x1[2], M1)
var x2i = Mod.call(A2[0]*_x2[0] + A2[1]*_x2[1] + A2[2]*_x2[2], M2)
_x1 = [x1i, _x1[0], _x1[1]] /* keep last three */
_x2 = [x2i, _x2[0], _x2[1]] /* keep last three */
return Mod.call(x1i - x2i, M1) + 1
}
nextFloat { nextInt / D }
}
var randomGen = MRG32k3a.new()
randomGen.seed(1234567)
for (i in 0..4) System.print(randomGen.nextInt)
var counts = List.filled(5, 0)
randomGen.seed(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:
1459213977 2827710106 4245671317 3877608661 2595287583 The counts for 100,000 repetitions are: 0 : 20002 1 : 20060 2 : 19948 3 : 20059 4 : 19931