all tasks

This commit is contained in:
Ingy döt Net 2013-04-11 01:07:29 -07:00
parent b83f433714
commit 68f8f3e56b
14735 changed files with 178959 additions and 0 deletions

View file

@ -0,0 +1,33 @@
A ''subtractive generator'' calculates a sequence of [[random number generator|random numbers]], where each number is congruent to the subtraction of two previous numbers from the sequence. The formula is
* <math>r_n = r_{(n - i)} - r_{(n - j)} \pmod m</math>
for some fixed values of <math>i</math>, <math>j</math> and <math>m</math>, all positive integers. Supposing that <math>i > j</math>, then the state of this generator is the list of the previous numbers from <math>r_{n - i}</math> to <math>r_{n - 1}</math>. Many states generate uniform random integers from <math>0</math> to <math>m - 1</math>, but some states are bad. A state, filled with zeros, generates only zeros. If <math>m</math> is even, then a state, filled with even numbers, generates only even numbers. More generally, if <math>f</math> is a factor of <math>m</math>, then a state, filled with multiples of <math>f</math>, generates only multiples of <math>f</math>.
All subtractive generators have some weaknesses. The formula correlates <math>r_n</math>, <math>r_{(n - i)}</math> and <math>r_{(n - j)}</math>; these three numbers are not independent, as true random numbers would be. Anyone who observes <math>i</math> consecutive numbers can predict the next numbers, so the generator is not cryptographically secure. The authors of ''Freeciv'' ([http://svn.gna.org/viewcvs/freeciv/trunk/utility/rand.c?view=markup 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 <math>r_{(n - i)} - r_{(n - j)}</math> is always between <math>-m</math> and <math>m</math>, so a program only needs to add <math>m</math> to negative numbers.
The choice of <math>i</math> and <math>j</math> affects the period of the generator. A popular choice is <math>i = 55</math> and <math>j = 24</math>, so the formula is
* <math>r_n = r_{(n - 55)} - r_{(n - 24)} \pmod m</math>
The subtractive generator from ''xpat2'' uses
* <math>r_n = r_{(n - 55)} - r_{(n - 24)} \pmod{10^9}</math>
The implementation is by J. Bentley and comes from program_tools/universal.c of [ftp://dimacs.rutgers.edu/pub/netflow/ the DIMACS (netflow) archive] at Rutgers University. It credits Knuth, [[wp:The Art of Computer Programming|''TAOCP'']], Volume 2, Section 3.2.2 (Algorithm A).
Bentley uses this clever algorithm to seed the generator.
# Start with a single <math>seed</math> in range <math>0</math> to <math>10^9 - 1</math>.
# Set <math>s_0 = seed</math> and <math>s_1 = 1</math>. The inclusion of <math>s_1 = 1</math> avoids some bad states (like all zeros, or all multiples of 10).
# Compute <math>s_2, s_3, ..., s_{54}</math> using the subtractive formula <math>s_n = s_{(n - 2)} - s_{(n - 1)} \pmod{10^9}</math>.
# Reorder these 55 values so <math>r_0 = s_{34}</math>, <math>r_1 = s_{13}</math>, <math>r_2 = s_{47}</math>, ..., <math>r_n = s_{(34 * (n + 1) \pmod{55})}</math>.
#* This is the same order as <math>s_0 = r_{54}</math>, <math>s_1 = r_{33}</math>, <math>s_2 = r_{12}</math>, ..., <math>s_n = r_{((34 * n) - 1 \pmod{55})}</math>.
#* This rearrangement exploits how 34 and 55 are relatively prime.
# Compute the next 165 values <math>r_{55}</math> to <math>r_{219}</math>. Store the last 55 values.
This generator yields the sequence <math>r_{220}</math>, <math>r_{221}</math>, <math>r_{222}</math> and so on. For example, if the seed is 292929, then the sequence begins with <math>r_{220} = 467478574</math>, <math>r_{221} = 512932792</math>, <math>r_{222} = 539453717</math>. By starting at <math>r_{220}</math>, 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 <math>r_n</math>. Any array or list would work; a [[ring buffer]] is ideal but not necessary.
Implement a subtractive generator that replicates the sequences from ''xpat2''.

View file

@ -0,0 +1,2 @@
---
note: Randomness

View file

@ -0,0 +1,11 @@
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;

View file

@ -0,0 +1,33 @@
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;

View file

@ -0,0 +1,14 @@
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;

View file

@ -0,0 +1,26 @@
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

View file

@ -0,0 +1,52 @@
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$;

View file

@ -0,0 +1,63 @@
// 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;
}

View file

@ -0,0 +1,43 @@
#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;
}

View file

@ -0,0 +1,23 @@
(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))))

View file

@ -0,0 +1,52 @@
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());
}

View file

@ -0,0 +1,100 @@
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)
}
}
}

View file

@ -0,0 +1,8 @@
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

View file

@ -0,0 +1,26 @@
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

View file

@ -0,0 +1,17 @@
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

View file

@ -0,0 +1,6 @@
$ jconsole
load'sg.ijs'
sg 2
467478574 512932792
sg 4
539453717 20349702 615542081 378707948

View file

@ -0,0 +1,9 @@
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]}

View file

@ -0,0 +1,34 @@
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

View file

@ -0,0 +1,2 @@
sgv=vector(55,i,random(10^9));sgi=1;
sg()=sgv[sgi=sgi%55+1]=(sgv[sgi]-sgv[(sgi+30)%55+1])%10^9

View file

@ -0,0 +1,31 @@
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;

View file

@ -0,0 +1,17 @@
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];

View file

@ -0,0 +1,26 @@
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);

View file

@ -0,0 +1,18 @@
(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) )

View file

@ -0,0 +1,2 @@
(subRandSeed 292929)
(do 7 (println (subRand)))

View file

@ -0,0 +1,42 @@
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()

View file

@ -0,0 +1,23 @@
/*REXX pgm uses a subtractive generator, creates a seq of random numbers*/
numeric digits 20; billion = 10**9; s.0 = 292929; s.1 = 1
cI = 55; cJ = 24; cP = 34
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.*/
/*──────────────────────────────────subroutines─────────────────────────*/
r: parse arg _; return r._
s: parse arg _; return s._
mod: procedure; arg a,b; return ((a // b) + b) // b

View file

@ -0,0 +1,37 @@
# 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 }

View file

@ -0,0 +1,35 @@
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]
}