Subtractive generator
You are encouraged to solve this task according to the task description, using any language you may know.
A subtractive generator calculates a sequence of random numbers, where each number is congruent to the subtraction of two previous numbers from the sequence. The formula is
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_n = r_{(n - i)} - r_{(n - j)} \pmod m}
for some fixed values of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle i} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle j} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle m} , all positive integers. Supposing that Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle i > j} , then the state of this generator is the list of the previous numbers from Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_{n - i}} to Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_{n - 1}} . Many states generate uniform random integers from Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 0} to Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle m - 1} , but some states are bad. A state, filled with zeros, generates only zeros. If Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle m} is even, then a state, filled with even numbers, generates only even numbers. More generally, if Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle f} is a factor of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle m} , then a state, filled with multiples of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle f} , generates only multiples of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle f} .
All subtractive generators have some weaknesses. The formula correlates Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_n} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_{(n - i)}} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_{(n - j)}} ; these three numbers are not independent, as true random numbers would be. Anyone who observes Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle i} consecutive numbers can predict the next numbers, so the generator is not cryptographically secure. The authors of Freeciv (utility/rand.c) and xpat2 (src/testit2.c) knew another problem: the low bits are less random than the high bits.
The subtractive generator has a better reputation than the linear congruential generator, perhaps because it holds more state. A subtractive generator might never multiply numbers: this helps where multiplication is slow. A subtractive generator might also avoid division: the value of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_{(n - i)} - r_{(n - j)}} is always between Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle -m} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle m} , so a program only needs to add Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle m} to negative numbers.
The choice of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle i} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle j} affects the period of the generator. A popular choice is Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle i = 55} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle j = 24} , so the formula is
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_n = r_{(n - 55)} - r_{(n - 24)} \pmod m}
The subtractive generator from xpat2 uses
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_n = r_{(n - 55)} - r_{(n - 24)} \pmod{10^9}}
The implementation is by J. Bentley and comes from program_tools/universal.c of the DIMACS (netflow) archive at Rutgers University. It credits Knuth, TAOCP, Volume 2, Section 3.2.2 (Algorithm A).
Bentley uses this clever algorithm to seed the generator.
- Start with a single Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle seed} in range Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 0} to Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 10^9 - 1} .
- Set Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle s_0 = seed} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle s_1 = 1} . The inclusion of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle s_1 = 1} avoids some bad states (like all zeros, or all multiples of 10).
- Compute Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle s_2, s_3, ..., s_{54}} using the subtractive formula Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle s_n = s_{(n - 2)} - s_{(n - 1)} \pmod{10^9}} .
- Reorder these 55 values so Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_0 = s_{34}}
, Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_1 = s_{13}}
, Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_2 = s_{47}}
, ..., Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_n = s_{(34 * (n + 1) \pmod{55})}}
.
- This is the same order as Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle s_0 = r_{54}} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle s_1 = r_{33}} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle s_2 = r_{12}} , ..., Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle s_n = r_{((34 * n) - 1 \pmod{55})}} .
- This rearrangement exploits how 34 and 55 are relatively prime.
- Compute the next 165 values Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_{55}} to Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_{219}} . Store the last 55 values.
This generator yields the sequence Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_{220}} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_{221}} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_{222}} and so on. For example, if the seed is 292929, then the sequence begins with Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_{220} = 467478574} , , . By starting at , this generator avoids a bias from the first numbers of the sequence. This generator must store the last 55 numbers of the sequence, so to compute the next . Any array or list would work; a ring buffer is ideal but not necessary.
Implement a subtractive generator that replicates the sequences from xpat2.
Bracmat
This is a translation of the C example.
<lang bracmat> 1000000000:?MOD; tbl$(state,55); 0:?si:?sj;
(subrand-seed=
i,j,p2
. 1:?p2
& mod$(!arg,!MOD):?(0$?state)
& 1:?i
& 21:?j
& whl
' ( !i:<55
& (!j:~<55&!j+-55:?j|)
& !p2:?(!j$?state)
& ( !arg+-1*!p2:?p2:<0
& !p2+!MOD:?p2
|
)
& !(!j$state):?arg
& !i+1:?i
& !j+21:?j
)
& 0:?s1:?i
& 24:?sj
& whl
' ( !i:<165
& subrand$
& !i+1:?i
));
(subrand=
x
. (!si:!sj&subrand-seed$0|)
& (!si:>0&!si+-1|54):?si
& (!sj:>0&!sj+-1|54):?sj
& ( !(!si$state)+-1*!(!sj$state):?x:<0
& !x+!MOD:?x
|
)
& !x:?(!si$?state));
(Main=
i
. subrand-seed$292929
& 0:?i
& whl
' ( !i:<10
& out$(subrand$)
& !i+1:?i
));
Main$; </lang>
Output: <lang> 467478574 512932792 539453717 20349702 615542081 378707948 933204586 824858649 506003769 380969305 </lang>
C
This is basically the same as the reference C code, only differs in that it's C89. <lang c>#include<stdio.h>
- define MOD 1000000000
int state[55], si = 0, sj = 0;
int subrand();
void subrand_seed(int p1) { int i, j, p2 = 1;
state[0] = p1 % MOD; for (i = 1, j = 21; i < 55; i++, j += 21) { if (j >= 55) j -= 55; state[j] = p2; if ((p2 = p1 - p2) < 0) p2 += MOD; p1 = state[j]; } si = 0; sj = 24; for (i = 0; i < 165; i++) subrand(); }
int subrand() { int x; if (si == sj) subrand_seed(0);
if (!si--) si = 54; if (!sj--) sj = 54; if ((x = state[si] - state[sj]) < 0) x += MOD;
return state[si] = x; }
int main() { subrand_seed(292929); int i; for (i = 0; i < 10; i++) printf("%d\n", subrand());
return 0; }</lang>
Common Lisp
<lang lisp>(defun sub-rand (state)
(let ((x (last state)) (y (last state 25))) ;; I take "circular buffer" very seriously (until some guru ;; points out it's utterly wrong thing to do) (setf (cdr x) state) (lambda () (setf x (cdr x)
y (cdr y) (car x) (mod (- (car x) (car y)) (expt 10 9))))))
- returns an RNG with Bentley seeding
(defun bentley-clever (seed)
(let ((s (list 1 seed)) f)
(dotimes (i 53)
(push (mod (- (cadr s) (car s)) (expt 10 9)) s))
(setf f (sub-rand
(loop for i from 1 to 55 collect (elt s (- 54 (mod (* 34 i) 55))))))
(dotimes (x 165) (funcall f)) f))
- test it (output same as everyone else's)
(let ((f (bentley-clever 292929)))
(dotimes (x 10) (format t "~a~%" (funcall f))))</lang>
D
<lang d>import std.stdio;
struct Subtractive {
enum MOD = 1_000_000_000; private int[55] state; private int si, sj;
this(in int p1) pure nothrow {
subrandSeed(p1);
}
void subrandSeed(int p1) pure nothrow {
int p2 = 1;
state[0] = p1 % MOD;
for (int i = 1, j = 21; i < 55; i++, j += 21) {
if (j >= 55)
j -= 55;
state[j] = p2;
if ((p2 = p1 - p2) < 0)
p2 += MOD;
p1 = state[j];
}
si = 0;
sj = 24;
foreach (i; 0 .. 165)
subrand();
}
int subrand() pure nothrow {
if (si == sj)
subrandSeed(0);
if (!si--)
si = 54;
if (!sj--)
sj = 54;
int x = state[si] - state[sj];
if (x < 0)
x += MOD;
return state[si] = x; }
}
void main() {
auto gen = Subtractive(292_929);
foreach (i; 0 .. 10)
writeln(gen.subrand());
}</lang> Output:
467478574 512932792 539453717 20349702 615542081 378707948 933204586 824858649 506003769 380969305
dc
<lang dc>[*
* (seed) lsx -- * Seeds the subtractive generator. * Uses register R to hold the state. *]sz
[
[* Fill ring buffer R[0] to R[54]. *]sz d 54:R SA [A = R[54] = seed]sz 1 d 33:R SB [B = R[33] = 1]sz 12 SC [C = index 12, into array R.]sz [55 -]SI [ [Loop until C is 54:]sz lA lB - d lC:R [R[C] = A - B]sz lB sA sB [Parallel let A = B and B = A - B]sz lC 34 + d 55 !>I d sC [C += 34 (mod 55)]sz 54 !=L ]d SL x [* Point R[55] and R[56] into ring buffer. *]sz 0 55:R [R[55] = index 0, of 55th last number.]sz 31 56:R [R[56] = index 31, of 24th last number.]sz [* Stir ring buffer. *]sz 165 [ [Loop 165 times:]sz 55;R;R 56;R;R - 55;R:R [Discard a random number.]sz 55;R 1 + d 55 !>I 55:R [R[55] += 1 (mod 55)]sz 56;R 1 + d 55 !>I 56:R [R[56] += 1 (mod 55)]sz 1 - d 0 <L ]d sL x LAsz LBsz LCsz LIsz LLsz
]ss
[*
* lrx -- (random number from 0 to 10^9 - 1) * Returns the next number from the subtractive generator. * Uses register R, seeded by lsx. *]sz
[
55;R;R 56;R;R - [R[R[55]] - R[R[56]] is next random number.]sz d 55;R:R [Put it in R[R[55]]. Also leave it on stack.]sz [55 -]SI 55;R 1 + d 55 !>I 55:R [R[55] += 1 (mod 55)]sz 56;R 1 + d 55 !>I 56:R [R[56] += 1 (mod 55)]sz [1000000000 +]sI 1000000000 % d 0 >I [Random number = it (mod 10^9)]sz LIsz
]sr
[* Seed with 292929 and print first three random numbers. *]sz
292929 lsx
lrx psz
lrx psz
lrx psz</lang>
This program prints 467478574, 512932792, 539453717.
This implementation never uses multiplication, but it does use modulus (remainder from division) to put each random number in range from 0 to 10^9 - 1.
Go
<lang go>package main
import (
"fmt" "os"
)
// A fairly close port of the Bentley code, but parameterized to better // conform to the algorithm description in the task, which didn't assume // constants for i, j, m, and seed. also parameterized here are k, // the reordering factor, and s, the number of intial numbers to discard, // as these are dependant on i. func newSG(i, j, k, s, m, seed int) func() int {
// check parameters for range and mutual consistency
assert(i > 0, "i must be > 0")
assert(j > 0, "j must be > 0")
assert(i > j, "i must be > j")
assert(k > 0, "k must be > 0")
p, q := i, k
if p < q {
p, q = q, p
}
for q > 0 {
p, q = q, p%q
}
assert(p == 1, "k, i must be relatively prime")
assert(s >= i, "s must be >= i")
assert(m > 0, "m must be > 0")
assert(seed >= 0, "seed must be >= 0")
// variables for closure f
arr := make([]int, i)
a := 0
b := j
// f is Bently RNG lprand
f := func() int {
if a == 0 {
a = i
}
a--
if b == 0 {
b = i
}
b--
t := arr[a] - arr[b]
if t < 0 {
t += m
}
arr[a] = t
return t
}
// Bentley seed algorithm sprand
last := seed
arr[0] = last
next := 1
for i0 := 1; i0 < i; i0++ {
ii := k * i0 % i
arr[ii] = next
next = last - next
if next < 0 {
next += m
}
last = arr[ii]
}
for i0 := i; i0 < s; i0++ {
f()
}
// return the fully initialized RNG
return f
}
func assert(p bool, m string) {
if !p {
fmt.Println(m)
os.Exit(1)
}
}
func main() {
// 1st test case included in program_tools/universal.c.
// (2nd test case fails. A single digit is missing, indicating a typo.)
ptTest(0, 1, []int{921674862, 250065336, 377506581})
// reproduce 3 values given in task description
skip := 220
sg := newSG(55, 24, 21, skip, 1e9, 292929)
for n := skip; n <= 222; n++ {
fmt.Printf("r(%d) = %d\n", n, sg())
}
}
func ptTest(nd, s int, rs []int) {
sg := newSG(55, 24, 21, 220+nd, 1e9, s)
for _, r := range rs {
a := sg()
if r != a {
fmt.Println("Fail")
os.Exit(1)
}
}
}</lang> Output:
r(220) = 467478574 r(221) = 512932792 r(222) = 539453717
Icon and Unicon
<lang Icon>procedure main()
every 1 to 10 do
write(rand_sub(292929))
end
procedure rand_sub(x) static ring,m
if /ring then {
m := 10^9
every (seed | ring) := list(55)
seed[1] := \x | ?(m-1)
seed[2] := 1
every seed[n := 3 to 55] := (seed[n-2]-seed[n-1])%m
every ring[(n := 0 to 54) + 1] := seed[1 + (34 * (n + 1)%55)]
every n := *ring to 219 do {
ring[1] -:= ring[-24]
ring[1] %= m
put(ring,get(ring))
}
}
ring[1] -:= ring[-24]
ring[1] %:= m
if ring[1] < 0 then ring[1] +:= m
put(ring,get(ring))
return ring[-1]
end</lang>
Output:
467478574 512932792 539453717 20349702 615542081 378707948 933204586 824858649 506003769 380969305
OCaml
<lang ocaml>let _mod = 1_000_000_000 let state = Array.create 55 0 let si = ref 0 let sj = ref 0
let rec subrand_seed _p1 =
let p1 = ref _p1 in let p2 = ref 1 in state.(0) <- !p1 mod _mod; let j = ref 21 in for i = 1 to pred 55 do if !j >= 55 then j := !j - 55; state.(!j) <- !p2; p2 := !p1 - !p2; if !p2 < 0 then p2 := !p2 + _mod; p1 := state.(!j); j := !j + 21; done; si := 0; sj := 24; for i = 0 to pred 165 do ignore (subrand()) done
and subrand() =
if !si = !sj then subrand_seed 0; decr si; if !si < 0 then si := 54; decr sj; if !sj < 0 then sj := 54; let x = state.(!si) - state.(!sj) in let x = if x < 0 then x + _mod else x in state.(!si) <- x; (x)
let () =
subrand_seed 292929; for i = 1 to 10 do Printf.printf "%d\n" (subrand()) done</lang>
Output:
$ ocaml sub_gen.ml 467478574 512932792 539453717 20349702 615542081 378707948 933204586 824858649 506003769 380969305
PARI/GP
<lang parigp>sgv=vector(55,i,random(10^9));sgi=1; sg()=sgv[sgi=sgi%55+1]=(sgv[sgi]-sgv[(sgi+30)%55+1])%10^9</lang>
Perl
<lang perl>use 5.10.0; use strict;
{ # bracket state data into a lexical scope my @state; my $mod = 1_000_000_000;
sub bentley_clever { my @s = ( shift() % $mod, 1); push @s, ($s[-2] - $s[-1]) % $mod while @s < 55; @state = map($s[(34 + 34 * $_) % 55], 0 .. 54); subrand() for (55 .. 219); }
sub subrand() { bentley_clever(0) unless @state; # just incase
my $x = (shift(@state) - $state[-24]) % $mod; push @state, $x; $x; } }
bentley_clever(292929);
say subrand() for (1 .. 10);</lang>output
467478574 512932792 539453717 20349702 615542081 ...
Perl 6
<lang perl6>sub bentley_clever($seed) {
constant $mod = 1_000_000_000; my @seeds = ($seed % $mod, 1, (* - *) % $mod ... *)[^55]; my @state = @seeds[ 34, (* + 34 ) % 55 ... 0 ];
sub subrand() {
push @state, (my $x = (@state.shift - @state[*-24]) % $mod);
$x;
}
subrand for 55 .. 219;
&subrand ... *;
}
my @sr := bentley_clever(292929); .say for @sr[^10];</lang> Here we just make the seeder return the random sequence as a lazy list.
Output:
467478574 512932792 539453717 20349702 615542081 378707948 933204586 824858649 506003769 380969305
PicoLisp
Using a circular list (as a true "ring" buffer). <lang PicoLisp>(setq
*Bentley (apply circ (need 55)) *Bentley2 (nth *Bentley 32) )
(de subRandSeed (S)
(let (N 1 P (nth *Bentley 55))
(set P S)
(do 54
(set (setq P (nth P 35)) N)
(when (lt0 (setq N (- S N)))
(inc 'N 1000000000) )
(setq S (car P)) ) )
(do 165 (subRand)) )
(de subRand ()
(when (lt0 (dec *Bentley (pop '*Bentley2)))
(inc *Bentley 1000000000) )
(pop '*Bentley) )</lang>
Test: <lang PicoLisp>(subRandSeed 292929) (do 7 (println (subRand)))</lang> Output:
467478574 512932792 539453717 20349702 615542081 378707948 933204586
Python
Uses collections.deque as a ring buffer
<lang python> import collections s= collections.deque(maxlen=55)
- Start with a single seed in range 0 to 10**9 - 1.
seed = 292929
- Set s0 = seed and s1 = 1.
- The inclusion of s1 = 1 avoids some bad states
- (like all zeros, or all multiples of 10).
s.append(seed) s.append(1)
- Compute s2,s3,...,s54 using the subtractive formula
- sn = s(n - 2) - s(n - 1)(mod 10**9).
for n in xrange(2, 55):
s.append((s[n-2] - s[n-1]) % 10**9)
- Reorder these 55 values so r0 = s34, r1 = s13, r2 = s47, ...,
- rn = s(34 * (n + 1)(mod 55)).
for i in [0,1,2,3,50,51,52,53,54]:
print "s[", i, "] = ", s[i]
r = collections.deque(maxlen=55) for n in xrange(55):
i = (34 * (n+1)) % 55 r.append(s[i])
- This is the same order as s0 = r54, s1 = r33, s2 = r12, ...,
- sn = r((34 * n) - 1(mod 55)).
- This rearrangement exploits how 34 and 55 are relatively prime.
- Compute the next 165 values r55 to r219. Store the last 55 values.
def getnextr():
"""get next random number""" r.append((r[0]-r[31])%10**9) return r[54]
- rn = r(n - 55) - r(n - 24)(mod 10**9) for n >= 55
for n in xrange(219 - 54):
getnextr()
- now fully initilised
- print first five numbers
for i in xrange(5):
print "result = ", getnextr()
</lang>
Ruby
This implementation aims for simplicity, not speed. SubRandom#rand pushes to and shifts from an array; this might be slower than a ring buffer. The seeding method must call rand 55 extra times (220 times instead of 165 times). The code also calls Ruby's modulus operator, which always returns a non-negative integer if the modulus is positive.
<lang ruby># SubRandom is a subtractive random number generator which generates
- the same sequences as Bentley's generator, as used in xpat2.
class SubRandom
# The original seed of this generator. attr_reader :seed
# Creates a SubRandom generator with the given _seed_.
# The _seed_ must be an integer from 0 to 999_999_999.
def initialize(seed = Kernel.rand(1_000_000_000))
(0..999_999_999).include? seed or
raise ArgumentError, "seed not in 0..999_999_999"
# @state = 55 elements.
ary = [seed, 1]
53.times { ary << ary[-2] - ary[-1] }
@state = []
34.step(1870, 34) {|i| @state << ary[i % 55] }
220.times { rand } # Discard first 220 elements of sequence.
@seed = seed # Save original seed. end
# Duplicates internal state so SubRandom#dup never shares state. def initialize_copy(orig) @state = @state.dup end
# Returns the next random integer, from 0 to 999_999_999. def rand @state << (@state[-55] - @state[-24]) % 1_000_000_000 @state.shift end
end
rng = SubRandom.new(292929) p (1..3).map { rng.rand }</lang>
[467478574, 512932792, 539453717]
Tcl
<lang tcl>package require Tcl 8.5 namespace eval subrand {
variable mod 1000000000 state [lrepeat 55 0] si 0 sj 0
proc seed p1 {
global subrand::mod subrand::state subrand::si subrand::sj set p2 1 lset state 0 [expr {$p1 % $mod}] for {set i 1; set j 21} {$i < 55} {incr i; incr j 21} { if {$j >= 55} {incr j -55} lset state $j $p2 if {[set p2 [expr {$p1 - $p2}]] < 0} {incr p2 $mod} set p1 [lindex $state $j] } set si 0 set sj 24 for {set i 0} {$i < 165} {incr i} { gen }
}
proc gen {} {
global subrand::mod subrand::state subrand::si subrand::sj if {$si == $sj} {seed 0} if {[incr si -1] < 0} {set si 54} if {[incr sj -1] < 0} {set sj 54} set x [expr {[lindex $state $si] - [lindex $state $sj]}] if {$x < 0} {incr x $mod} lset state $si $x return $x
}
}
subrand::seed 292929 for {set i 0} {$i < 10} {incr i} {
puts [subrand::gen]
}</lang>