(* random.sml * * COPYRIGHT (c) 1993 by AT&T Bell Laboratories. See COPYRIGHT file for details. * * This package implements a random number generator using a subtract-with-borrow * (SWB) generator as described in Marsaglia and Zaman, "A New Class of Random Number * Generators," Ann. Applied Prob. 1(3), 1991, pp. 462-480. * * The SWB generator is a 31-bit generator with lags 48 and 8. It has period * (2^1487 - 2^247)/105 or about 10^445. In general, these generators are * excellent. However, they act locally like a lagged Fibonacci generator * and thus have troubles with the birthday test. Thus, we combine this SWB * generator with the linear congruential generator (48271*a)mod(2^31-1). * * Although the interface is fairly abstract, the implementation uses * 31-bit ML words. At some point, it might be good to use 32-bit words. *) structure Random : RANDOM = struct structure A = Array structure LW = LargeWord structure W8A = Word8Array structure W8V = Word8Vector structure P = Pack32Big val << = Word31.<< val >> = Word31.>> val & = Word31.andb val ++ = Word31.orb val xorb = Word31.xorb infix << >> & ++ val nbits = 31 (* bits per word *) val maxWord : Word31.word = 0wx7FFFFFFF (* largest word *) val bit30 : Word31.word = 0wx40000000 val lo30 : Word31.word = 0wx3FFFFFFF val N = 48 val lag = 8 val offset = N-lag fun error (f,msg) = LibBase.failure {module="Random",func=f, msg=msg} val two2neg30 = 1.0/((real 0x8000)*(real 0x8000)) (* 2^~30 *) fun minus(x,y,false) = (x - y, y > x) | minus(x,y,true) = (x - y - 0w1, y >= x) datatype rand = RND of { vals : Word31.word A.array,(* seed array *) borrow : bool ref, (* last borrow *) congx : Word31.word ref, (* congruential seed *) index : int ref (* index of next available value in vals *) } (* We represent state as a string, starting with an initial * word acting as an magic cookie (with bit 0 determining the * value of borrow), followed by a word containing index and a word * containing congx, followed by the seed array. *) val numWords = 3 + N val magic : LW.word = 0wx72646e64 fun toString (RND{vals, borrow, congx, index}) = let val arr = W8A.array (4*numWords, 0w0) val word0 = if !borrow then LW.orb (magic, 0w1) else magic fun fill (src,dst) = if src = N then () else ( P.update (arr, dst, Word31.toLargeWord (A.sub (vals, src))); fill (src+1,dst+1) ) in P.update (arr, 0, word0); P.update (arr, 1, LW.fromInt (!index)); P.update (arr, 2, Word31.toLargeWord (!congx)); fill (0,3); Byte.bytesToString (W8A.extract (arr, 0, NONE)) end fun fromString s = let val bytes = Byte.stringToBytes s val _ = if W8V.length bytes = 4 * numWords then () else error ("fromString","invalid state string") val word0 = P.subVec (bytes, 0) val _ = if LW.andb(word0, 0wxFFFFFFFE) = magic then () else error ("fromString","invalid state string") fun subVec i = P.subVec (bytes, i) val borrow = ref (LW.andb(word0,0w1) = 0w1) val index = ref (LW.toInt (subVec 1)) val congx = ref (Word31.fromLargeWord (subVec 2)) val arr = A.array (N, 0w0 : Word31.word) fun fill (src,dst) = if dst = N then () else ( A.update (arr, dst, Word31.fromLargeWord (subVec src)); fill (src+1,dst+1) ) in fill (3, 0); RND{vals = arr, index = index, congx = congx, borrow = borrow} end (* linear congruential generator: * multiplication by 48271 mod (2^31 - 1) *) val a : Word31.word = 0w48271 val m : Word31.word = 0w2147483647 val q = m div a val r = m mod a fun lcg seed = let val left = a * (seed mod q) val right = r * (seed div q) in if left > right then left - right else (m - right) + left end (* Fill seed array using subtract-with-borrow generator: * x[n] = x[n-lag] - x[n-N] - borrow * Sets index to 1 and returns 0th value. *) fun fill (RND{vals,index,congx,borrow}) = let fun update (ix,iy,b) = let val (z,b') = minus(A.sub(vals,ix), A.sub(vals,iy),b) in A.update(vals,iy,z); b' end fun fillup (i,b) = if i = lag then b else fillup(i+1, update(i+offset,i,b)) fun fillup' (i,b) = if i = N then b else fillup'(i+1, update(i-lag,i,b)) in borrow := fillup' (lag, fillup (0,!borrow)); index := 1; A.sub(vals,0) end (* Create initial seed array and state of generator. * Fills the seed array one bit at a time by taking the leading * bit of the xor of a shift register and a congruential sequence. * The congruential generator is (c*48271) mod (2^31 - 1). * The shift register generator is c(I + L18)(I + R13). * The same congruential generator continues to be used as a * mixing generator with the SWB generator. *) fun rand (congy, shrgx) = let fun mki (i,c,s) = let val c' = lcg c val s' = xorb(s, s << 0w18) val s'' = xorb(s', s' >> 0w13) val i' = (lo30 & (i >> 0w1)) ++ (bit30 & xorb(c',s'')) in (i',c',s'') end fun iterate (0, v) = v | iterate (n, v) = iterate(n-1, mki v) fun mkseed (congx,shrgx) = iterate (nbits, (0w0,congx,shrgx)) fun genseed (0,seeds,congx,_) = (seeds,congx) | genseed (n,seeds,congx,shrgx) = let val (seed,congx',shrgx') = mkseed (congx,shrgx) in genseed(n-1,seed::seeds,congx',shrgx') end val congx = ((Word31.fromInt congy & maxWord) << 0w1)+0w1 val (seeds,congx) = genseed(N,[],congx, Word31.fromInt shrgx) in RND{vals = A.fromList seeds, index = ref 0, congx = ref congx, borrow = ref false} end (* Get next random number. The tweak function combines * the number from the SWB generator with a number from * the linear congruential generator. *) fun randWord (r as RND{vals, index,congx,...}) = let val idx = !index fun tweak i = let val c = lcg (!congx) in congx := c; xorb(i, c) end in if idx = N then tweak(fill r) else tweak(A.sub(vals,idx)) before index := idx+1 end fun randInt state = Word31.toIntX(randWord state) fun randNat state = Word31.toIntX(randWord state & lo30) fun randReal state = (real(randNat state) + real(randNat state) * two2neg30) * two2neg30 fun randRange (i,j) = if j < i then error ("randRange", "hi < lo") else let val R = two2neg30*real(j - i + 1) in fn s => i + trunc(R*real(randNat s)) end handle _ => let val ri = real i val R = (real j)-ri+1.0 in fn s => trunc(ri + R*(randReal s)) end end; (* Random *)
Click to toggle
does not end with </html> tag
does not end with </body> tag
The output has ended thus: in fn s => trunc(ri + R*(randReal s)) end end; (* Random *)