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
)




(
mod

m
)



{\displaystyle r_{n}=r_{(n-i)}-r_{(n-j)}{\pmod {m}}}


for some fixed values of 



i


{\displaystyle i}

, 



j


{\displaystyle j}

 and 



m


{\displaystyle m}

, all positive integers. Supposing that 



i
>
j


{\displaystyle i>j}

, then the state of this generator is the list of the previous numbers from 




r

n
−
i




{\displaystyle r_{n-i}}

 to 




r

n
−
1




{\displaystyle r_{n-1}}

. Many states generate uniform random integers from 



0


{\displaystyle 0}

 to 



m
−
1


{\displaystyle m-1}

, but some states are bad. A state, filled with zeros, generates only zeros. If 



m


{\displaystyle m}

 is even, then a state, filled with even numbers, generates only even numbers. More generally, if 



f


{\displaystyle f}

 is a factor of 



m


{\displaystyle m}

, then a state, filled with multiples of 



f


{\displaystyle f}

, generates only multiples of 



f


{\displaystyle f}

.

All subtractive generators have some weaknesses. The formula correlates 




r

n




{\displaystyle r_{n}}

, 




r

(
n
−
i
)




{\displaystyle r_{(n-i)}}

 and 




r

(
n
−
j
)




{\displaystyle r_{(n-j)}}

; these three numbers are not independent, as true random numbers would be. Anyone who observes 



i


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




r

(
n
−
i
)


−

r

(
n
−
j
)




{\displaystyle r_{(n-i)}-r_{(n-j)}}

 is always between 



−
m


{\displaystyle -m}

 and 



m


{\displaystyle m}

, so a program only needs to add 



m


{\displaystyle m}

 to negative numbers.

The choice of 



i


{\displaystyle i}

 and 



j


{\displaystyle j}

 affects the period of the generator. A popular choice is 



i
=
55


{\displaystyle i=55}

 and 



j
=
24


{\displaystyle j=24}

, so the formula is

 




r

n


=

r

(
n
−
55
)


−

r

(
n
−
24
)




(
mod

m
)



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


The subtractive generator from xpat2 uses

 




r

n


=

r

(
n
−
55
)


−

r

(
n
−
24
)




(
mod


10

9


)



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



s
e
e
d


{\displaystyle seed}

 in range 



0


{\displaystyle 0}

 to 




10

9


−
1


{\displaystyle 10^{9}-1}

.
 Set 




s

0


=
s
e
e
d


{\displaystyle s_{0}=seed}

 and 




s

1


=
1


{\displaystyle s_{1}=1}

. The inclusion of 




s

1


=
1


{\displaystyle s_{1}=1}

 avoids some bad states (like all zeros, or all multiples of 10).
 Compute 




s

2


,

s

3


,
.
.
.
,

s

54




{\displaystyle s_{2},s_{3},...,s_{54}}

 using the subtractive formula 




s

n


=

s

(
n
−
2
)


−

s

(
n
−
1
)




(
mod


10

9


)



{\displaystyle s_{n}=s_{(n-2)}-s_{(n-1)}{\pmod {10^{9}}}}

.
 Reorder these 55 values so 




r

0


=

s

34




{\displaystyle r_{0}=s_{34}}

, 




r

1


=

s

13




{\displaystyle r_{1}=s_{13}}

, 




r

2


=

s

47




{\displaystyle r_{2}=s_{47}}

, ..., 




r

n


=

s

(
34
∗
(
n
+
1
)


(
mod

55
)

)




{\displaystyle r_{n}=s_{(34*(n+1){\pmod {55}})}}

.
 This is the same order as 




s

0


=

r

54




{\displaystyle s_{0}=r_{54}}

, 




s

1


=

r

33




{\displaystyle s_{1}=r_{33}}

, 




s

2


=

r

12




{\displaystyle s_{2}=r_{12}}

, ..., 




s

n


=

r

(
(
34
∗
n
)
−
1


(
mod

55
)

)




{\displaystyle 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




{\displaystyle r_{55}}

 to 




r

219




{\displaystyle r_{219}}

. Store the last 55 values.
This generator yields the sequence 




r

220




{\displaystyle r_{220}}

, 




r

221




{\displaystyle r_{221}}

, 




r

222




{\displaystyle r_{222}}

 and so on. For example, if the seed is 292929, then the sequence begins with 




r

220


=
467478574


{\displaystyle r_{220}=467478574}

, 




r

221


=
512932792


{\displaystyle r_{221}=512932792}

, 




r

222


=
539453717


{\displaystyle r_{222}=539453717}

. By starting at 




r

220




{\displaystyle 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




{\displaystyle 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.




i,j,m,d,seed = 55,24,10^9,34,292929     # parameters
s = Array{Int32}(undef,i); r = similar(s)
s[1:2] = [seed,1]                       # table initialization
for n = 3:i; (s[n] = s[n-2]-s[n-1]) < 0 && (s[n] += m) end
t = 1; for u=1:i; (global t+=d)>i && (t-=i); r[u]=s[t] end # permutation, r = s[(d*(1:i) .% i).+1]
 
u,v,n = i,i-j,i-1
while (n += 1) > 0
    (global u += 1) > i && (u = 1)      # circular indexing: u,v = ((n,n-j) .% i).+1
    (global v += 1) > i && (v = 1)
    (r[u] -= r[v]) < 0 && (r[u] += m)   # table update
    n < 220 && continue                 # 165 silent values
    print((n,r[u]))                     # show (index,value) of next pseudorandom number
    x = readline(stdin)                 # wait until the ENTER key is pressed
    length(x) > 0 && break              # any other key before ENTER => exit
end