diff options
Diffstat (limited to 'lib/std/rand.myr')
| -rw-r--r-- | lib/std/rand.myr | 185 |
1 files changed, 185 insertions, 0 deletions
diff --git a/lib/std/rand.myr b/lib/std/rand.myr new file mode 100644 index 0000000..2e8c73e --- /dev/null +++ b/lib/std/rand.myr @@ -0,0 +1,185 @@ +use "die.use" +use "types.use" +use "alloc.use" +/* + Translated from C by Ori Bernstein + */ + +/* + A C-program for MT19937, with initialization improved 2002/1/26. + Coded by Takuji Nishimura and Makoto Matsumoto. + + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + Any feedback is very welcome. + http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) + */ + +pkg std = + type rng + + const mksrng : (seed : uint32 -> rng#) + const delrng : (rng : rng# -> void) + generic rand : (rng : rng#, lo : @a::(numeric,integral), hi : @a::(numeric,integral) -> @a::(numeric,integral)) + generic randN : (rng : rng# -> @a::(numeric,integral)) + const randbytes : (rng : rng#, buf : byte[:] -> size) + const rand32 : (rng : rng# -> uint32) +;; + +type rng = struct + state : uint32[624] + i : uint32 +;; + +/* allocates and initializes a random number generator */ +const mksrng = {seed + var rng + + rng = alloc() + init(rng, seed) + -> rng +} + +const delrng = {rng + free(rng) +} + +/* initializes a random number generator from the seed `seed`. */ +const init = {rng, seed + var i + + for i = 0; i < 624; i++ + rng.state[i] = seed + seed = 1812433253 * (seed ^ (seed >> 30)) + i + 1 + ;; + rng.i = i +} + +/* + Generates a random integer from `rng` in the range [lo, hi), + returning the value. The range [lo, hi) must be positive, + nonempty, and the difference between hi and lo must be + less then 2^(type_bits - 1) +*/ +generic rand = {rng, lo, hi -> @a::(integral,numeric) + var span, lim + var maxrand + var val + + assert(hi - lo > 0, "rand.myr: range for random values must be >= 1") + + span = hi - lo + maxrand = (1 << (8*sizeof(@a))) - 1 /* max for signed value */ + if maxrand < 0 /* signed */ + maxrand = (1 << (8*sizeof(@a)-1)) - 1 /* max for signed value */ + ;; + + lim = (maxrand/span)*span + val = (randN(rng) & maxrand) + while val > lim + val = (randN(rng) & maxrand) + ;; + -> val % span + lo +} + +/* + Generates a random integer of any size from the + random number generator `rng`. The returned value + may be negative, if the type is signed. +*/ +generic randN = {rng -> @a::(integral,numeric) + var i, val + + val = 0 + for i = 0; i < sizeof(@a)/4; i++ + val <<= 8*sizeof(@a) + val |= rand32(rng) castto(@a::(integral,numeric)) + ;; + -> val +} + +/* + generates a 32 bit unsigned random number + from the random number generator `rng`. +*/ +const rand32 = {rng + var x + + if rng.i == 624 + next(rng) + ;; + x = rng.state[rng.i] + rng.i++ + + x ^= x >> 11 + x ^= (x << 7) & 0x9D2C5680 + x ^= (x << 15) & 0xEFC60000 + -> x ^ (x >> 18) +} + +const randbytes = {rng, buf + var i, n, r + + n = 0 + for i = 0; i < buf.len/4; i++ + r = rand32(rng) + buf[n++] = (r >> 0 & 0xff) castto(byte) + buf[n++] = (r >> 8 & 0xff) castto(byte) + buf[n++] = (r >> 16 & 0xff) castto(byte) + buf[n++] = (r >> 32 & 0xff) castto(byte) + ;; + r = rand32(rng) + for ; n != buf.len; n++ + buf[n++] = (r & 0xff) castto(byte) + r >>= 8 + ;; + -> n + +} + +/* updates random number generator state when we tick over. */ +const next = {rng + var k + var y + + for k = 0; k < 227; k++ + y = (rng.state[k] & 0x80000000) | (rng.state[k + 1] & 0x7FFFFFFF) + rng.state[k] = rng.state[k + 397] ^ (y >> 1) ^ ((y & 1) * 0x9908B0DF) + ;; + for ; k < 623; k++ + y = (rng.state[k] & 0x80000000) | (rng.state[k + 1] & 0x7FFFFFFF) + rng.state[k] = rng.state[k - 227] ^ (y >> 1) ^ ((y & 1) * 0x9908B0DF); + ;; + y = (rng.state[623] & 0x80000000) | (rng.state[0] & 0x7FFFFFFF) + rng.state[623] = rng.state[396] ^ (y >> 1) ^ ((y & 1) * 0x9908B0DF); + rng.i = 0 +} |