ref: 458b9914044d8ddcdaea19e875a76493494b60f0
parent: 93123de228f52f9294dcba6748d8e8e9a8446335
author: Ori Bernstein <ori@eigenstate.org>
date: Sat Sep 10 09:28:44 EDT 2016
Clean up the RNG. Now it's less complex/shitty.
--- a/lib/std/rand.myr
+++ b/lib/std/rand.myr
@@ -1,53 +1,12 @@
+use "alloc"
+use "assert"
use "die"
use "extremum"
-use "assert"
-use "types"
-use "alloc"
+use "mk"
use "now"
+use "putint"
+use "types"
-/*
- 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
@@ -56,153 +15,91 @@
generic rand : (lo : @a::(numeric,integral), hi : @a::(numeric,integral) -> @a::(numeric,integral))
generic randnum : (-> @a::(numeric,integral))
- const randbytes : (buf : byte[:] -> size)
+ const randbytes : (buf : byte[:] -> void)
generic rngrand : (rng : rng#, lo : @a::(numeric,integral), hi : @a::(numeric,integral) -> @a::(numeric,integral))
generic rngrandnum : (rng : rng# -> @a::(numeric,integral))
- const rngrandbytes : (rng : rng#, buf : byte[:] -> size)
+ const rngrandbytes : (rng : rng#, buf : byte[:] -> void)
;;
type rng = struct
- state : uint32[624]
- i : uint32
+ s0 : uint64
+ s1 : uint64
;;
-var _rng : rng#
+var _rng
+generic rand = {lo, hi; -> rngrand(&_rng, lo, hi)}+generic randnum = {; -> rngrandnum(&_rng)}+const randbytes = {buf; -> rngrandbytes(&_rng, buf)}+
const __init__ = {- _rng = mksrng((now() : uint32))
+ _rng.s0 = (now() : uint64)
+ _rng.s1 = (now() : uint64)
}
-/* allocates and initializes a random number generator */
-const mksrng = {seed- var rng
-
- rng = alloc()
- init(rng, seed)
- -> rng
+const mksrng = {seed : uint32 -> rng#+ -> std.mk([
+ .s0=(seed & 0xffff : uint64),
+ .s1=(seed&0xffff0000>>16 : uint64)
+ ])
}
-const freerng = {rng- free(rng)
+const freerng = {r+ std.free(r)
}
-/* initializes a random number generator from the seed `seed`. */
-const init = {rng, seed- for var i = 0; i < 624; i++
- rng.state[i] = seed
- seed = 1812433253 * (seed ^ (seed >> 30)) + i + 1
- ;;
- rng.i = 624
-}
-
-generic rand = {lo, hi- -> rngrand(_rng, lo, hi)
-}
-
-generic randnum = {- -> rngrandnum(_rng)
-}
-
-const randbytes = {buf- -> rngrandbytes(_rng, buf)
-}
-
-/*
- 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 rngrand = {rng, lo, hi -> @a::(integral,numeric)- var span, lim
- var maxrand
- var val
-
+generic rngrand = {rng, lo, hi+ var span, lim, val, max
+
span = abs(hi - lo)
- assert(span > 0, "rand.myr: range for random values must be >= 1, got hi={}, lo={}\n", hi, lo)-
- maxrand = ~0
- if maxrand < 0 /* if we're signed the maximum value is different */
- maxrand = (1 << (8*sizeof(@a)-1)) - 1 /* max for signed value */
+ max = ~0
+ /* if ~0 is negative, we have a signed value with a different max */
+ if max < 0
+ max = (1 << (8*sizeof(@a)-1)) - 1
;;
-
- lim = (maxrand/span)*span
- val = (rngrandnum(rng) & maxrand)
- while val > lim
- val = (rngrandnum(rng) & maxrand)
+
+ lim = (max/span)*span
+ val = (rngrandnum(rng) & max)
+ while val > lim
+ val = (rngrandnum(rng) & max)
;;
-> 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 rngrandnum = {rng -> @a::(integral,numeric)- var val
-
- val = 0
- for var i = 0; i < sizeof(@a)/4; i++
- val <<= 8*sizeof(@a)
- val |= (rand32(rng) : @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 rngrandbytes = {rng, buf- var n, r
+ var n, r : uint64
n = 0
- for var i = 0; i < buf.len/4; i++
- r = rand32(rng)
- buf[n++] = (r >> 0 & 0xff : byte)
- buf[n++] = (r >> 8 & 0xff : byte)
- buf[n++] = (r >> 16 & 0xff : byte)
- buf[n++] = (r >> 32 & 0xff : byte)
+ for var i = 0; i + 8 < buf.len/8; i++
+ r = rngrandnum(rng)
+ putle64(buf[n:n+8], r)
+ n += 8
;;
- r = rand32(rng)
+ r = rngrandnum(rng)
for ; n != buf.len; n++
buf[n] = (r & 0xff : 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
+/*
+Generate a number using the xoroshiro algorithm, as
+designed by David Blackman and Sebastiano Vigna.
+
+See http://xoroshiro.di.unimi.it/ for details.
+*/
+generic rngrandnum = {rng -> @a::(numeric,integral)+ var s0, s1, r
+
+ s0 = rng.s0
+ s1 = rng.s1
+
+ r = s0 + s1
+ s1 ^= s0
+
+ rng.s0 = (s0 << 55 | s0 >> 9) ^ s1 ^ (s1 << 14)
+ rng.s1 = (s1 << 36 | s1 >> 28)
+ -> (r : @a::(numeric,integral))
}