Subtractive generator
From Rosetta Code
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
- rn = r(n − i) − r(n − j)(mod m)
for some fixed values of i, j and m, all positive integers. Supposing that i > j, then the state of this generator is the list of the previous numbers from rn − i to rn − 1. Many states generate uniform random integers from 0 to m − 1, but some states are bad. A state, filled with zeros, generates only zeros. If m is even, then a state, filled with even numbers, generates only even numbers. More generally, if f is a factor of m, then a state, filled with multiples of f, generates only multiples of f.
All subtractive generators have some weaknesses. The formula correlates rn, r(n − i) and r(n − j); these three numbers are not independent, as true random numbers would be. Anyone who observes 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 r(n − i) − r(n − j) is always between − m and m, so a program only needs to add m to negative numbers.
The choice of i and j affects the period of the generator. A popular choice is i = 55 and j = 24, so the formula is
- rn = r(n − 55) − r(n − 24)(mod m)
The subtractive generator from xpat2 uses
- rn = r(n − 55) − r(n − 24)(mod 109)
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 seed in range 0 to 109 − 1.
- Set s0 = seed and s1 = 1. The inclusion of s1 = 1 avoids some bad states (like all zeros, or all multiples of 10).
- Compute s2,s3,...,s54 using the subtractive formula sn = s(n − 2) − s(n − 1)(mod 109).
- Reorder these 55 values so r0 = s34, r1 = s13, r2 = s47, ..., rn = s(34 * (n + 1)(mod 55)).
- 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.
This generator yields the sequence r220, r221, r222 and so on. For example, if the seed is 292929, then the sequence begins with r220 = 467478574, r221 = 512932792, r222 = 539453717. By starting at r220, 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 rn. Any array or list would work; a ring buffer is ideal but not necessary.
Implement a subtractive generator that replicates the sequences from xpat2.
Contents |
[edit] Ada
subtractive_generator.ads:
package Subtractive_Generator is
type State is private;
procedure Initialize (Generator : in out State; Seed : Natural);
procedure Next (Generator : in out State; N : out Natural);
private
type Number_Array is array (Natural range <>) of Natural;
type State is record
R : Number_Array (0 .. 54);
Last : Natural;
end record;
end Subtractive_Generator;
subtractive_generator.adb:
package body Subtractive_Generator is
procedure Initialize (Generator : in out State; Seed : Natural) is
S : Number_Array (0 .. 1);
I : Natural := 0;
J : Natural := 1;
begin
S (0) := Seed;
S (1) := 1;
Generator.R (54) := S (0);
Generator.R (33) := S (1);
for N in 2 .. Generator.R'Last loop
S (I) := (S (I) - S (J)) mod 10 ** 9;
Generator.R ((34 * N - 1) mod 55) := S (I);
I := (I + 1) mod 2;
J := (J + 1) mod 2;
end loop;
Generator.Last := 54;
for I in 1 .. 165 loop
Subtractive_Generator.Next (Generator => Generator, N => J);
end loop;
end Initialize;
procedure Next (Generator : in out State; N : out Natural) is
begin
Generator.Last := (Generator.Last + 1) mod 55;
Generator.R (Generator.Last) :=
(Generator.R (Generator.Last)
- Generator.R ((Generator.Last - 24) mod 55)) mod 10 ** 9;
N := Generator.R (Generator.Last);
end Next;
end Subtractive_Generator;
Example main.adb:
with Ada.Text_IO;
with Subtractive_Generator;
procedure Main is
Random : Subtractive_Generator.State;
N : Natural;
begin
Subtractive_Generator.Initialize (Generator => Random,
Seed => 292929);
for I in 220 .. 222 loop
Subtractive_Generator.Next (Generator => Random, N => N);
Ada.Text_IO.Put_Line (Integer'Image (I) & ":" & Integer'Image (N));
end loop;
end Main;
Output:
220: 467478574 221: 512932792 222: 539453717
[edit] Bracmat
This is a translation of the C example.
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$;
Output:
467478574 512932792 539453717 20349702 615542081 378707948 933204586 824858649 506003769 380969305
[edit] C
This is basically the same as the reference C code, only differs in that it's C89.
#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;
}
[edit] C++
// written for clarity not efficiency.
#include <iostream>
using std::cout;
using std::endl;
#include <boost/array.hpp>
#include <boost/circular_buffer.hpp>
class Subtractive_generator {
private:
static const int param_i = 55;
static const int param_j = 24;
static const int initial_load = 219;
static const int mod = 1e9;
boost::circular_buffer<int> r;
public:
Subtractive_generator(int seed);
int next();
int operator()(){return next();}
};
Subtractive_generator::Subtractive_generator(int seed)
:r(param_i)
{
boost::array<int, param_i> s;
s[0] = seed;
s[1] = 1;
for(int n = 2; n < param_i; ++n){
int t = s[n-2]-s[n-1];
if (t < 0 ) t+= mod;
s[n] = t;
}
for(int n = 0; n < param_i; ++n){
int i = (34 * (n+1)) % param_i;
r.push_back(s[i]);
}
for(int n = param_i; n <= initial_load; ++n) next();
}
int Subtractive_generator::next()
{
int t = r[0]-r[31];
if (t < 0) t += mod;
r.push_back(t);
return r[param_i-1];
}
int main()
{
Subtractive_generator rg(292929);
cout << "result = " << rg() << endl;
cout << "result = " << rg() << endl;
cout << "result = " << rg() << endl;
cout << "result = " << rg() << endl;
cout << "result = " << rg() << endl;
cout << "result = " << rg() << endl;
cout << "result = " << rg() << endl;
return 0;
}
Output:
result = 467478574 result = 512932792 result = 539453717 result = 20349702 result = 615542081 result = 378707948 result = 933204586
[edit] Common 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))))
[edit] 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());
}
Output:
467478574 512932792 539453717 20349702 615542081 378707948 933204586 824858649 506003769 380969305
[edit] 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
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.
[edit] 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)
}
}
}
Output:
r(220) = 467478574 r(221) = 512932792 r(222) = 539453717
[edit] Icon and Unicon
procedure main()Output:
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
467478574 512932792 539453717 20349702 615542081 378707948 933204586 824858649 506003769 380969305
[edit] J
sg.ijs
Loops are hidden in a generalized power conjunction ^: . f^:n y evaluates f n times, as in f(f(f(...f(y)))...) . Yes! f^:(-1) IS the inverse of f . When known.
came_from_locale_sg_=: coname''
cocurrent'sg' NB. install the state of rng sg into locale sg
SEED=: 292929
'I J M first_Bentley_number B2'=: 55 24 1e9 34 165
SG=: 1 : 'M&|@:-/@:(m&{)'
r=: (I|(first_Bentley_number*>:i.I)) { (, _2 _1 SG)^:(I-2) 1,~SEED
sg=: 3 : 0
t=. (, (-I,J)SG)^:y r
r=: y }. t
t {.~ -y
)
discard=. sg B2
cocurrent came_from_locale NB. return to previous locale
sg=: sg_sg_ NB. make a local name for sg in locale sg
Use:
$ jconsole
load'sg.ijs'
sg 2
467478574 512932792
sg 4
539453717 20349702 615542081 378707948
[edit] 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
Output:
$ ocaml sub_gen.ml 467478574 512932792 539453717 20349702 615542081 378707948 933204586 824858649 506003769 380969305
[edit] PARI/GP
sgv=vector(55,i,random(10^9));sgi=1;
sg()=sgv[sgi=sgi%55+1]=(sgv[sgi]-sgv[(sgi+30)%55+1])%10^9
[edit] Perl
use 5.10.0;output
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);
467478574512932792 539453717 20349702 615542081
...
[edit] Perl 6
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];
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
[edit] PicoLisp
Using a circular list (as a true "ring" buffer).
(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) )
Test:
(subRandSeed 292929)
(do 7 (println (subRand)))
Output:
467478574 512932792 539453717 20349702 615542081 378707948 933204586
[edit] Python
Uses collections.deque as a ring buffer
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)).
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()
[edit] 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.
# 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 }
[467478574, 512932792, 539453717]
[edit] 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]
}
