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.
Ada
subtractive_generator.ads: <lang Ada>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;</lang>
subtractive_generator.adb: <lang Ada>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;</lang>
Example main.adb: <lang Ada>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;</lang>
Output:
220: 467478574 221: 512932792 222: 539453717
AutoHotkey
<lang AutoHotkey>r := InitR(292929)
Loop, 10 Out .= (A_Index + 219) ":`t" GetRand(r) "`n"
MsgBox, % Out
GetRand(r) { i := Mod(r["j"], 55) , r[i] := Mod(r[i] - r[Mod(i + 31, 55)], r["m"]) , r["j"] += 1 return, (r[i] < 0 ? r[i] + r["m"] : r[i]) }
InitR(Seed) { r := {"j": 0, "m": 10 ** 9}, s := {0: Seed, 1: 1} Loop, 53 s[A_Index + 1] := Mod(s[A_Index - 1] - s[A_Index], r["m"]) Loop, 55 r[A_Index - 1] := s[Mod(34 * A_Index, 55)] Loop, 165 i := Mod(A_Index + 54, 55) , r[i] := Mod(r[i] - r[Mod(A_Index + 30, 55)], r["m"]) return, r }</lang> Output:
220: 467478574 221: 512932792 222: 539453717 223: 20349702 224: 615542081 225: 378707948 226: 933204586 227: 824858649 228: 506003769 229: 380969305
BBC BASIC
<lang bbcbasic> dummy% = FNsubrand(292929)
FOR i% = 1 TO 10 PRINT FNsubrand(0) NEXT END DEF FNsubrand(s%) PRIVATE r%(), p% : DIM r%(54) IF s% = 0 THEN p% = (p% + 1) MOD 55 r%(p%) = r%(p%) - r%((p% + 31) MOD 55) IF r%(p%) < 0 r%(p%) += 10^9 = r%(p%) ENDIF LOCAL i% r%(54) = s% : r%(33) = 1 p% = 12 FOR i% = 2 TO 54 r%(p%) = r%((p%+42) MOD 55) - r%((p%+21) MOD 55) IF r%(p%) < 0 r%(p%) += 10^9 p% = (p% + 34) MOD 55 NEXT FOR i% = 55 TO 219 IF FNsubrand(0) NEXT = 0</lang>
Output:
467478574 512932792 539453717 20349702 615542081 378707948 933204586 824858649 506003769 380969305
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:
467478574 512932792 539453717 20349702 615542081 378707948 933204586 824858649 506003769 380969305
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>
C++
<lang cpp> // 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;
} </lang>
Output:
result = 467478574 result = 512932792 result = 539453717 result = 20349702 result = 615542081 result = 378707948 result = 933204586
C#
<lang csharp> public class SubtractiveGenerator {
public static int MAX = 1000000000; private int[] state; private int pos;
private int mod(int n) { return ((n % MAX) + MAX) % MAX; }
public SubtractiveGenerator(int seed) { state = new int[55];
int[] temp = new int[55]; temp[0] = mod(seed); temp[1] = 1; for(int i = 2; i < 55; ++i) temp[i] = mod(temp[i - 2] - temp[i - 1]);
for(int i = 0; i < 55; ++i) state[i] = temp[(34 * (i + 1)) % 55];
pos = 54; for(int i = 55; i < 220; ++i) next(); }
public int next() { int temp = mod(state[(pos + 1) % 55] - state[(pos + 32) % 55]); pos = (pos + 1) % 55; state[pos] = temp; return temp; }
static void Main(string[] args) { SubtractiveGenerator gen = new SubtractiveGenerator(292929); for(int i = 220; i < 230; ++i) Console.WriteLine(i.ToString() + ": " + gen.next().ToString()); }
} </lang>
Output:
220: 467478574 221: 512932792 222: 539453717 223: 20349702 224: 615542081 225: 378707948 226: 933204586 227: 824858649 228: 506003769 229: 380969305
Clojure
<lang clojure>(defn xpat2-with-seed
"produces an xpat2 function initialized from seed" [seed] (let [e9 1000000000 fs (fn i j [j (mod (- i j) e9)]) s (->> [seed 1] (iterate fs) (map first) (take 55) vec) rinit (map #(-> % inc (* 34) (mod 55) s) (range 55)) r-atom (atom [54 (int-array rinit)]) update (fn nprev r (let [n (-> nprev inc (mod 55)) rx #(get r (-> n (- %) (mod 55))) rn (-> (rx 55) (- (rx 24)) (mod e9)) _ (aset-int r n rn)] [n r])) xpat2 #(let [[n r] (swap! r-atom update)] (get r n)) _ (dotimes [_ 165] (xpat2))] xpat2))
(def xpat2 (xpat2-with-seed 292929))
(println (xpat2) (xpat2) (xpat2)) ; prints: 467478574 512932792 539453717 </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
Haskell
<lang haskell>subtractgen seed = drop 220 out where out = mmod $ r ++ zipWith (-) out (drop 31 out) where r = take 55 $ shuffle $ cycle $ take 55 s shuffle x = head xx:shuffle xx where xx = drop 34 x s = mmod $ seed:1:zipWith (-) s (tail s) mmod = map (`mod` 10^9)
main = mapM_ print $ take 10 $ subtractgen 292929</lang>
- Output:
467478574 512932792 539453717 20349702 615542081 378707948 933204586 824858649 506003769 380969305
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
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.
<lang J>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 </lang>
Use: <lang sh>$ jconsole
load'sg.ijs' sg 2
467478574 512932792
sg 4
539453717 20349702 615542081 378707948
</lang>
Mathematica
<lang Mathematica>initialize[n_] :=
Module[{buffer}, buffer = Join[Nest[Flatten@{#, Mod[Subtract @@ #-2 ;;, 10^9]} &, {n, 1}, 53][[1 + Mod[34 Range@54, 55]]], {n}]; Nest[nextValue, buffer, 165]] nextValue[buffer_] := Flatten@{Rest@buffer, Mod[Subtract @@ buffer[[{1, 32}]], 10^9]}</lang>
buffer = initialize[292929]; Do[Print@Last[buffer = nextValue[buffer]], {10}] 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
PL/I
<lang PL/I> subtractive_generator: procedure options (main);
declare (r, s) (0:54) fixed binary (31); declare (i, n, seed) fixed binary (31);
/* Bentley's initialization */ seed = 292929; s(0) = seed; s(1) = 1;
/* Compute s2,s3,...,s54 using the subtractive formula sn = s(n-2) - s(n-1)(mod 10**9). */ do n = 2 to hbound(s,1); s(n) = mod ( s(n-2) - s(n-1), 1000000000); end;
/* Rearrange initial values. */ do n = 0 to hbound(r,1); r(n) = s( mod(34*(n+1), 55)); end;
do n = 55 to 219; i = mod (n, 55); r(i) = mod ( r(mod(n-55, 55)) - r(mod(n-24, 55)), 1000000000); end;
do n = 220 to 235; i = mod(n, 55); r(i) = mod ( r(mod(n-55, 55)) - r(mod(n-24, 55)), 1000000000); put skip list (r(i)); end;
end subtractive_generator; </lang>
Required 3 results: 467478574 512932792 539453717 Subsequent values: 20349702 615542081 378707948 933204586 824858649 506003769 380969305 442823364 994162810 261423281 139610325 80746560 563900213
Python
Python: With explanation
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)).
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>
Python: As a class within a module
Python 2 and 3 compatable. <lang python>import collections
_ten2nine = 10**9
class Subtractive_generator():
def __init__(self, seed=292929): self.r = collections.deque(maxlen=55) s = collections.deque(maxlen=55) s.extend([seed, 1]) s.extend((s[n-2] - s[n-1]) % _ten2nine for n in range(2, 55)) self.r.extend(s[(34 * (n+1)) % 55] for n in range(55)) for n in range(219 - 54): self() def __call__(self): r = self.r r.append((r[0] - r[31]) % _ten2nine) return r[54]
if __name__ == '__main__':
srand = Subtractive_generator() print([srand() for i in range(5)])</lang>
- Output:
[467478574, 512932792, 539453717, 20349702, 615542081]
Racket
<lang Racket>#lang racket (define (make-initial-state a-list max-i)
(for/fold ((state a-list)) ((i (in-range (length a-list) max-i))) (append state (list (- (list-ref state (- i 2)) (list-ref state (- i 1))))))) ;from the seed and 1 creates the initial state
(define (shuffle a-list)
(for/list ((i (in-range (length a-list)))) (list-ref a-list (modulo (* 34 (add1 i)) 55)))) ;shuffles the state
(define (advance-state state (times 1))
(cond ((= 0 times) state) (else (advance-state (cdr (append state (list (modulo (- (list-ref state 0) (list-ref state 31)) (expt 10 9))))) (sub1 times))))) ;takes a state and the times it must be advanced, and returns the new state
(define (create-substractive-generator s0)
(define s1 1) (define first-state (make-initial-state (list s0 s1) 55)) (define shuffled-state (shuffle first-state)) (define last-state (advance-state shuffled-state 165)) (lambda ((m (expt 10 9))) (define new-state (advance-state last-state)) (set! last-state new-state) (modulo (car (reverse last-state)) m))) ;the lambda is a function with an optional argument ;that returns a new random number each time it's called
(define rand (create-substractive-generator 292929)) (build-list 3 (lambda (_) (rand))) ;returns a list made from the 3 wanted numbers</lang>
REXX
This REXX program is essentially a copy of the PL/I version. <lang rexx>/*REXX pgm uses a subtractive generator, creates a sequence of random #s*/ numeric digits 20; s.0=292929; s.1=1; billion=10**9 cI=55; cJ=24; cP=34; billion=1e9 /*same*/
do i=2 to cI-1 s.i=mod(s(i-2) - s(i-1), billion) end /*i*/ do j=0 to cI-1 r.j=s(mod(cP*(j+1), cI)) end /*j*/
m=219
do k=cI to m; x=k//cI r.x=mod(r(mod(k-cI, cI)) - r(mod(k-cJ, cI)), billion) end /*m*/
t=235
do n=m+1 to t; y=n//cI r.y=mod(r(mod(n-cI, cI)) - r(mod(n-cJ, cI)), billion) say right(r.y, 40) end /*n*/
exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────one─liner subroutines───────────────*/ mod: procedure; parse arg a,b; return ((a // b) + b) // b r: parse arg _; return r._ s: parse arg _; return s._</lang> output when using the default input:
467478574 512932792 539453717 20349702 615542081 378707948 933204586 824858649 506003769 380969305 442823364 994162810 261423281 139610325 80746560 563900213
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]
Seed7
<lang seed7>$ include "seed7_05.s7i";
const integer: MOD is 1000000000;
const type: subtractiveGenerator is new struct
var array integer: state is [0 .. 54] times 0; var integer: si is 0; var integer: sj is 24; end struct;
const func integer: subrand (inout subtractiveGenerator: generator) is forward;
const func subtractiveGenerator: subrandSeed (in var integer: p1) is func
result var subtractiveGenerator: generator is subtractiveGenerator.value; local var integer: p2 is 1; var integer: i is 0; var integer: j is 21; begin generator.state[0] := p1 mod MOD; for i range 1 to 54 do generator.state[j] := p2; p2 := (p1 - p2) mod MOD; p1 := generator.state[j]; j := (j + 21) mod 55; end for; for i range 1 to 165 do ignore(subrand(generator)); end for; end func;
const func integer: subrand (inout subtractiveGenerator: generator) is func
result var integer: subrand is 0; begin if generator.si = generator.sj then generator := subrandSeed(0); end if; generator.si := pred(generator.si) mod 55; generator.sj := pred(generator.sj) mod 55; subrand := (generator.state[generator.si] - generator.state[generator.sj]) mod MOD; generator.state[generator.si] := subrand; end func;
const proc: main is func
local var subtractiveGenerator: gen is subrandSeed(292929); var integer: i is 0; begin for i range 1 to 10 do writeln(subrand(gen)); end for; end func;</lang>
- Output:
467478574 512932792 539453717 20349702 615542081 378707948 933204586 824858649 506003769 380969305
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>
uBasic/4tH
<lang>Push 292929 : Gosub 100 : d = Pop()
For i = 1 To 10
Push 0 : Gosub 100 Print Pop()
Next
End
100 s = Pop()
If s = 0 Then p = (p + 1) % 55 @(p) = @(p) - @((p + 31) % 55) If @(p) < 0 Then @(p) = @(p) + 1000000000 Endif Push (@(p)) : Return Endif
@(54) = s : @(33) = 1 p = 12
For i = 2 To 54 @(p) = @((p + 42) % 55) - @((p + 21) % 55) If @(p) < 0 Then @(p) = @(p) + 1000000000 Endif p = (p + 34) % 55 Next
For i = 55 To 219 Push 0 : Gosub 100 : d = Pop() Next
Push 0 : Return</lang>
Output:
467478574 512932792 539453717 20349702 615542081 378707948 933204586 824858649 506003769 380969305 0 OK, 0:864
zkl
<lang zkl>fcn rand_sub(x){
var ring=L(),m=(1e9).toInt(); mod:='wrap(n){ if(n<0) n+m else n }; if(not ring){ seed:=L( (if(vm.numArgs) x else m-1), 1); foreach n in ([2 .. 54]){ seed.append((seed[n-2]-seed[n-1]):mod(_)) } foreach n in (55){ ring.append(seed[(34*(n+1))%55]) } do(220-ring.len()){ self.fcn() } // 165 } ring.append((ring.pop(0)-ring[-24]):mod(_)); return(ring[-1]);
}</lang> <lang zkl>do(4){ println(rand_sub(292929)) } //seed ignored after first call</lang>
- Output:
467478574 512932792 539453717 20349702