summaryrefslogtreecommitdiff
path: root/lib/std/rand.myr
diff options
context:
space:
mode:
authorOri Bernstein <ori@eigenstate.org>2016-09-10 13:28:44 -0700
committerOri Bernstein <ori@eigenstate.org>2016-09-10 13:30:41 -0700
commit458b9914044d8ddcdaea19e875a76493494b60f0 (patch)
tree09e107d4db28d7fc4c8e60f8954dc873ebe6427d /lib/std/rand.myr
parent93123de228f52f9294dcba6748d8e8e9a8446335 (diff)
downloadmc-458b9914044d8ddcdaea19e875a76493494b60f0.tar.gz
Clean up the RNG.
Now it's less complex/shitty.
Diffstat (limited to 'lib/std/rand.myr')
-rw-r--r--lib/std/rand.myr221
1 files changed, 59 insertions, 162 deletions
diff --git a/lib/std/rand.myr b/lib/std/rand.myr
index fce743c..c2a03fd 100644
--- a/lib/std/rand.myr
+++ b/lib/std/rand.myr
@@ -1,52 +1,11 @@
+use "alloc"
+use "assert"
use "die"
use "extremum"
-use "assert"
-use "types"
-use "alloc"
+use "mk"
use "now"
-
-/*
- 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)
- */
+use "putint"
+use "types"
pkg std =
type rng
@@ -56,153 +15,91 @@ pkg std =
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#
-
-const __init__ = {
- _rng = mksrng((now() : uint32))
-}
+var _rng
-/* allocates and initializes a random number generator */
-const mksrng = {seed
- var rng
+generic rand = {lo, hi; -> rngrand(&_rng, lo, hi)}
+generic randnum = {; -> rngrandnum(&_rng)}
+const randbytes = {buf; -> rngrandbytes(&_rng, buf)}
- rng = alloc()
- init(rng, seed)
- -> rng
-}
-
-const freerng = {rng
- free(rng)
-}
-
-/* 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)
+const __init__ = {
+ _rng.s0 = (now() : uint64)
+ _rng.s1 = (now() : uint64)
}
-generic randnum = {
- -> rngrandnum(_rng)
+const mksrng = {seed : uint32 -> rng#
+ -> std.mk([
+ .s0=(seed & 0xffff : uint64),
+ .s1=(seed&0xffff0000>>16 : uint64)
+ ])
}
-const randbytes = {buf
- -> rngrandbytes(_rng, buf)
+const freerng = {r
+ std.free(r)
}
-/*
- 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))
}