Home My Page Projects Code Snippets Project Openings SML/NJ
Summary Activity Forums Tracker Lists Tasks Docs Surveys News SCM Files

SCM Repository

[smlnj] Annotation of /sml/trunk/src/smlnj-lib/Util/random.sml
ViewVC logotype

Annotation of /sml/trunk/src/smlnj-lib/Util/random.sml

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1350 - (view) (download)

1 : monnier 2 (* random.sml
2 :     *
3 :     * COPYRIGHT (c) 1993 by AT&T Bell Laboratories. See COPYRIGHT file for details.
4 :     *
5 :     * This package implements a random number generator using a subtract-with-borrow
6 :     * (SWB) generator as described in Marsaglia and Zaman, "A New Class of Random Number
7 :     * Generators," Ann. Applied Prob. 1(3), 1991, pp. 462-480.
8 :     *
9 :     * The SWB generator is a 31-bit generator with lags 48 and 8. It has period
10 :     * (2^1487 - 2^247)/105 or about 10^445. In general, these generators are
11 :     * excellent. However, they act locally like a lagged Fibonacci generator
12 :     * and thus have troubles with the birthday test. Thus, we combine this SWB
13 :     * generator with the linear congruential generator (48271*a)mod(2^31-1).
14 :     *
15 :     * Although the interface is fairly abstract, the implementation uses
16 :     * 31-bit ML words. At some point, it might be good to use 32-bit words.
17 :     *)
18 :    
19 :     structure Random : RANDOM =
20 :     struct
21 :     structure A = Array
22 :     structure LW = LargeWord
23 :     structure W8A = Word8Array
24 :     structure W8V = Word8Vector
25 :     structure P = Pack32Big
26 :    
27 :     val << = Word31.<<
28 :     val >> = Word31.>>
29 :     val & = Word31.andb
30 :     val ++ = Word31.orb
31 :     val xorb = Word31.xorb
32 :     infix << >> & ++
33 :    
34 :     val nbits = 31 (* bits per word *)
35 :     val maxWord : Word31.word = 0wx7FFFFFFF (* largest word *)
36 :     val bit30 : Word31.word = 0wx40000000
37 :     val lo30 : Word31.word = 0wx3FFFFFFF
38 :    
39 :     val N = 48
40 :     val lag = 8
41 :     val offset = N-lag
42 :    
43 :     fun error (f,msg) = LibBase.failure {module="Random",func=f, msg=msg}
44 :    
45 :     val two2neg30 = 1.0/((real 0x8000)*(real 0x8000)) (* 2^~30 *)
46 :    
47 :     fun minus(x,y,false) = (x - y, y > x)
48 :     | minus(x,y,true) = (x - y - 0w1, y >= x)
49 :    
50 :     datatype rand = RND of {
51 :     vals : Word31.word A.array,(* seed array *)
52 :     borrow : bool ref, (* last borrow *)
53 :     congx : Word31.word ref, (* congruential seed *)
54 :     index : int ref (* index of next available value in vals *)
55 :     }
56 :    
57 :     (* We represent state as a string, starting with an initial
58 :     * word acting as an magic cookie (with bit 0 determining the
59 :     * value of borrow), followed by a word containing index and a word
60 :     * containing congx, followed by the seed array.
61 :     *)
62 :     val numWords = 3 + N
63 :     val magic : LW.word = 0wx72646e64
64 :     fun toString (RND{vals, borrow, congx, index}) = let
65 :     val arr = W8A.array (4*numWords, 0w0)
66 :     val word0 = if !borrow then LW.orb (magic, 0w1) else magic
67 :     fun fill (src,dst) =
68 :     if src = N then ()
69 :     else (
70 :     P.update (arr, dst, Word31.toLargeWord (A.sub (vals, src)));
71 :     fill (src+1,dst+1)
72 :     )
73 :     in
74 :     P.update (arr, 0, word0);
75 :     P.update (arr, 1, LW.fromInt (!index));
76 :     P.update (arr, 2, Word31.toLargeWord (!congx));
77 :     fill (0,3);
78 : mblume 1350 Byte.bytesToString (W8A.vector arr)
79 : monnier 2 end
80 :    
81 :     fun fromString s = let
82 :     val bytes = Byte.stringToBytes s
83 :     val _ = if W8V.length bytes = 4 * numWords then ()
84 :     else error ("fromString","invalid state string")
85 :     val word0 = P.subVec (bytes, 0)
86 :     val _ = if LW.andb(word0, 0wxFFFFFFFE) = magic then ()
87 :     else error ("fromString","invalid state string")
88 :     fun subVec i = P.subVec (bytes, i)
89 :     val borrow = ref (LW.andb(word0,0w1) = 0w1)
90 :     val index = ref (LW.toInt (subVec 1))
91 :     val congx = ref (Word31.fromLargeWord (subVec 2))
92 :     val arr = A.array (N, 0w0 : Word31.word)
93 :     fun fill (src,dst) =
94 :     if dst = N then ()
95 :     else (
96 :     A.update (arr, dst, Word31.fromLargeWord (subVec src));
97 :     fill (src+1,dst+1)
98 :     )
99 :     in
100 :     fill (3, 0);
101 :     RND{vals = arr,
102 :     index = index,
103 :     congx = congx,
104 :     borrow = borrow}
105 :     end
106 :    
107 :     (* linear congruential generator:
108 :     * multiplication by 48271 mod (2^31 - 1)
109 :     *)
110 :     val a : Word31.word = 0w48271
111 :     val m : Word31.word = 0w2147483647
112 :     val q = m div a
113 :     val r = m mod a
114 :     fun lcg seed = let
115 :     val left = a * (seed mod q)
116 :     val right = r * (seed div q)
117 :     in
118 :     if left > right then left - right
119 :     else (m - right) + left
120 :     end
121 :    
122 :     (* Fill seed array using subtract-with-borrow generator:
123 :     * x[n] = x[n-lag] - x[n-N] - borrow
124 :     * Sets index to 1 and returns 0th value.
125 :     *)
126 :     fun fill (RND{vals,index,congx,borrow}) = let
127 :     fun update (ix,iy,b) = let
128 :     val (z,b') = minus(A.sub(vals,ix), A.sub(vals,iy),b)
129 :     in
130 :     A.update(vals,iy,z); b'
131 :     end
132 :     fun fillup (i,b) =
133 :     if i = lag then b
134 :     else fillup(i+1, update(i+offset,i,b))
135 :     fun fillup' (i,b) =
136 :     if i = N then b
137 :     else fillup'(i+1, update(i-lag,i,b))
138 :     in
139 :     borrow := fillup' (lag, fillup (0,!borrow));
140 :     index := 1;
141 :     A.sub(vals,0)
142 :     end
143 :    
144 :     (* Create initial seed array and state of generator.
145 :     * Fills the seed array one bit at a time by taking the leading
146 :     * bit of the xor of a shift register and a congruential sequence.
147 :     * The congruential generator is (c*48271) mod (2^31 - 1).
148 :     * The shift register generator is c(I + L18)(I + R13).
149 :     * The same congruential generator continues to be used as a
150 :     * mixing generator with the SWB generator.
151 :     *)
152 :     fun rand (congy, shrgx) = let
153 :     fun mki (i,c,s) = let
154 :     val c' = lcg c
155 :     val s' = xorb(s, s << 0w18)
156 :     val s'' = xorb(s', s' >> 0w13)
157 :     val i' = (lo30 & (i >> 0w1)) ++ (bit30 & xorb(c',s''))
158 :     in (i',c',s'') end
159 :     fun iterate (0, v) = v
160 :     | iterate (n, v) = iterate(n-1, mki v)
161 :     fun mkseed (congx,shrgx) = iterate (nbits, (0w0,congx,shrgx))
162 :     fun genseed (0,seeds,congx,_) = (seeds,congx)
163 :     | genseed (n,seeds,congx,shrgx) = let
164 :     val (seed,congx',shrgx') = mkseed (congx,shrgx)
165 :     in genseed(n-1,seed::seeds,congx',shrgx') end
166 :     val congx = ((Word31.fromInt congy & maxWord) << 0w1)+0w1
167 :     val (seeds,congx) = genseed(N,[],congx, Word31.fromInt shrgx)
168 :     in
169 :     RND{vals = A.fromList seeds,
170 :     index = ref 0,
171 :     congx = ref congx,
172 :     borrow = ref false}
173 :     end
174 :    
175 :     (* Get next random number. The tweak function combines
176 :     * the number from the SWB generator with a number from
177 :     * the linear congruential generator.
178 :     *)
179 :     fun randWord (r as RND{vals, index,congx,...}) = let
180 :     val idx = !index
181 :     fun tweak i = let
182 :     val c = lcg (!congx)
183 :     in
184 :     congx := c;
185 :     xorb(i, c)
186 :     end
187 :     in
188 :     if idx = N then tweak(fill r)
189 :     else tweak(A.sub(vals,idx)) before index := idx+1
190 :     end
191 :    
192 :     fun randInt state = Word31.toIntX(randWord state)
193 :     fun randNat state = Word31.toIntX(randWord state & lo30)
194 :     fun randReal state =
195 :     (real(randNat state) + real(randNat state) * two2neg30) * two2neg30
196 :    
197 :     fun randRange (i,j) =
198 :     if j < i
199 :     then error ("randRange", "hi < lo")
200 :     else let
201 :     val R = two2neg30*real(j - i + 1)
202 :     in
203 :     fn s => i + trunc(R*real(randNat s))
204 :     end handle _ => let
205 :     val ri = real i
206 :     val R = (real j)-ri+1.0
207 :     in
208 :     fn s => trunc(ri + R*(randReal s))
209 :     end
210 :    
211 :     end; (* Random *)
212 :    

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0