Subtractive generator: Difference between revisions

Content added Content deleted

Inline

Revision as of 19:47, 4 January 2012

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

$r_{n}=r_{(n-i)}-r_{(n-j)}{\pmod {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 $r_{n-i}$ to $r_{n-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 $r_{n}$ , $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

$r_{n}=r_{(n-55)}-r_{(n-24)}{\pmod {m}}$

The subtractive generator from xpat2 uses

$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 $seed$ in range $0$ to $10^{9}-1$ .
Set $s_{0}=seed$ and $s_{1}=1$ . The inclusion of $s_{1}=1$ avoids some bad states (like all zeros, or all multiples of 10).
Compute $s_{2},s_{3},...,s_{54}$ using the subtractive formula $s_{n}=s_{(n-2)}-s_{(n-1)}{\pmod {10^{9}}}$ .
Reorder these 55 values so $r_{0}=s_{34}$ $r_0 = s_{34}$ , $r_{1}=s_{13}$ $r_1 = s_{13}$ , $r_{2}=s_{47}$ $r_2 = s_{47}$ , ..., $r_{n}=s_{(34*(n+1){\pmod {55}})}$ $r_n = s_{(34 * (n + 1) \pmod{55})}$ .
- This is the same order as $s_{0}=r_{54}$ , $s_{1}=r_{33}$ , $s_{2}=r_{12}$ , ..., $s_{n}=r_{((34*n)-1{\pmod {55}})}$ .
- This rearrangement exploits how 34 and 55 are relatively prime.
Compute the next 165 values $r_{55}$ to $r_{219}$ . Store the last 55 values.

This generator yields the sequence $r_{220}$ , $r_{221}$ , $r_{222}$ and so on. For example, if the seed is 292929, then the sequence begins with $r_{220}=467478574$ , $r_{221}=512932792$ , $r_{222}=539453717$ . By starting at $r_{220}$ , 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 $r_{n}$ . 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=

. (!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=

. subrand-seed$292929

 & 0:?i
 &   whl
   ' ( !i:<10
     & out$(subrand$)
     & !i+1:?i
     )
 & ok);

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

Translation of: C

<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:

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:

OCaml

Translation of: C

<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

Perl 6

Translation of: Perl

Works with: niecza

Works with: rakudo version nom

<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:

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:

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)).

print "s = ", s 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):

   print 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

Translation of: C

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