## random.pkg
# Compiled by:
# src/lib/std/standard.lib# 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, p 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 Mythryl words. At some point, it might be good to use 32-bit words. XXX BUGGO FIXME
#
######################################################################
# This is totally obsolete technology; we should
# implement something modern like
# http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/index.html
######################################################################
### "Anyone who considers arithmetic methods
### of producing random digits is, of course,
### in a state of sin."
### -- Johnny von Neuman
package random
: (weak) Random # Random is from src/lib/src/random.api{
package a = rw_vector; # rw_vector is from src/lib/std/src/rw-vector.pkg package lw = large_unt; # large_unt is from src/lib/std/large-unt.pkg package w8a = rw_vector_of_one_byte_unts; # rw_vector_of_one_byte_unts is from src/lib/std/src/rw-vector-of-one-byte-unts.pkg package w8v = vector_of_one_byte_unts; # vector_of_one_byte_unts is from src/lib/std/src/vector-of-one-byte-unts.pkg package p = pack_big_endian_unt1; # pack_big_endian_unt1 is from src/lib/std/src/pack-big-endian-unt1.pkg my (<<) = (tagged_unt::(<<));
my (>>) = (tagged_unt::(>>));
my (&) = tagged_unt::bitwise_and;
my (|) = tagged_unt::bitwise_or;
bitwise_xor = tagged_unt::bitwise_xor;
nbits = 31; # Bits per word
my max_word: tagged_unt::Unt = 0ux7FFFFFFF; # largest word
my bit30: tagged_unt::Unt = 0ux40000000;
my lo30: tagged_unt::Unt = 0ux3FFFFFFF;
nnn = 48;
lag = 8;
offset = nnn-lag;
fun error (f, msg)
=
lib_base::failure { module=>"Random", fn=>f, msg };
two2neg30 = 1.0 / ((float(0x8000))*(float(0x8000))); # 2^-30
Random_Number_Generator
=
RANDOM_NUMBER_GENERATOR {
vals: a::Rw_Vector( tagged_unt::Unt ), # Seed rw_vector
borrow: Ref( Bool ), # Last borrow
congx: Ref( tagged_unt::Unt ), # Congruential seed
index: Ref( Int ) # 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 rw_vector.
num_words = 3 + nnn;
my magic: lw::Unt
= 0ux72646e64;
fun to_string (RANDOM_NUMBER_GENERATOR { vals, borrow, congx, index } )
=
{ arr = w8a::make_rw_vector (4*num_words, 0u0);
word0 = if *borrow lw::bitwise_or (magic, 0u1);
else magic; fi;
fun fill (src, dst)
=
if (src != nnn)
p::set (arr, dst, tagged_unt::to_large_unt (a::get (vals, src)));
fill (src+1, dst+1);
fi;
p::set (arr, 0, word0);
p::set (arr, 1, lw::from_int *index);
p::set (arr, 2, tagged_unt::to_large_unt *congx);
fill (0, 3);
byte::bytes_to_string (w8a::to_vector arr);
};
fun from_string s
=
{ bytes = byte::string_to_bytes s;
if (w8v::length bytes != 4 * num_words)
error ("from_string", "invalid state string");
fi;
word0 = p::get_vec (bytes, 0);
if (lw::bitwise_and (word0, 0uxFFFFFFFE) != magic)
error ("from_string", "invalid state string");
fi;
fun get_vec i
=
p::get_vec (bytes, i);
borrow = REF (lw::bitwise_and (word0, 0u1) == 0u1);
index = REF (lw::to_int (get_vec 1));
congx = REF (tagged_unt::from_large_unt (get_vec 2));
arr = a::make_rw_vector (nnn, 0u0: tagged_unt::Unt);
fun fill (src, dst)
=
if (dst != nnn)
a::set (arr, dst, tagged_unt::from_large_unt (get_vec src));
fill (src+1, dst+1);
fi;
fill (3, 0);
RANDOM_NUMBER_GENERATOR
{ vals => arr,
index,
congx,
borrow
};
};
# linear congruential generator:
# multiplication by 48271 mod (2^31 - 1)
my a: tagged_unt::Unt = 0u48271;
my m: tagged_unt::Unt = 0u2147483647;
q = m / a;
r = m % a;
fun lcg seed
=
{ left = a * (seed % q);
right = r * (seed / q);
if (left > right) left - right;
else (m - right) + left; fi;
};
# Fill seed rw_vector using subtract-with-borrow generator:
# x[n] = x[n-lag] - x[n-nnn] - borrow
# Sets index to 1 and returns 0th value.
fun fill (RANDOM_NUMBER_GENERATOR { vals, index, congx, borrow } )
=
{
fun minus (x, y, FALSE) => (x - y, y > x);
minus (x, y, TRUE ) => (x - y - 0u1, y >= x);
end;
fun update (ix, iy, b)
=
{ my (z, b')
=
minus ( a::get (vals, ix),
a::get (vals, iy),
b
);
a::set (vals, iy, z); b';
};
fun fillup (i, b)
=
if (i == lag) b;
else fillup (i+1, update (i+offset, i, b));
fi;
fun fillup' (i, b)
=
if (i == nnn) b;
else fillup'(i+1, update (i-lag, i, b));
fi;
borrow := fillup' (lag, fillup (0,*borrow));
index := 1;
a::get (vals, 0);
};
# Create initial seed rw_vector and state of generator.
# Fills the seed rw_vector 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 make_random_number_generator (congy, shrgx)
=
{ fun mki (i, c, s)
=
{ c' = lcg c;
s' = bitwise_xor (s, s << 0u18);
s'' = bitwise_xor (s', s' >> 0u13);
i' = (lo30 & (i >> 0u1)) | (bit30 & bitwise_xor (c', s''));
(i', c', s'');
};
fun iterate (0, v) => v;
iterate (n, v) => iterate (n - 1, mki v);
end;
fun make_seed (congx, shrgx)
=
iterate (nbits, (0u0, congx, shrgx));
fun genseed (0, seeds, congx, _)
=>
(seeds, congx);
genseed (n, seeds, congx, shrgx)
=>
{ my (seed, congx', shrgx') = make_seed (congx, shrgx);
genseed (n - 1, seed ! seeds, congx', shrgx');
};
end;
congx = ((tagged_unt::from_int congy & max_word) << 0u1)+0u1;
my (seeds, congx)
=
genseed (nnn,[], congx, tagged_unt::from_int shrgx);
RANDOM_NUMBER_GENERATOR { vals => a::from_list seeds,
index => REF 0,
congx => REF congx,
borrow => REF FALSE };
};
# Get next random number. The tweak function combines
# the number from the SWB generator with a number from
# the linear congruential generator.
#
fun rand_word (r as RANDOM_NUMBER_GENERATOR { vals, index, congx, ... } )
=
{ idx = *index;
fun tweak i
=
{ c = lcg *congx;
congx := c;
bitwise_xor (i, c);
};
if (idx == nnn) tweak (fill r);
else tweak (a::get (vals, idx))
then
index := idx+1;
fi;
};
fun int state = tagged_unt::to_int_x (rand_word state);
fun nonnegative_int state = tagged_unt::to_int_x (rand_word state & lo30);
fun float01 state
=
(float(nonnegative_int state) + float(nonnegative_int state) * two2neg30) * two2neg30;
fun range (i, j)
=
if (j < i )
error ("random_range", "hi < lo");
else
{ rrr = two2neg30 * float(j - i + 1);
\\ s = i + trunc (rrr * float(nonnegative_int s));
}
except
_ = { ri = float(i);
rrr = (float(j))-ri+1.0;
\\ s = trunc (ri + rrr*(float01 s));
};
fi;
fun bool state
=
(int state) > 0;
}; # random