Commit dce7c00b authored by John van Groningen's avatar John van Groningen

64 bit port and also faster

parent a243f14b
......@@ -9,6 +9,8 @@
// generator" by M. Matsumoto and T. Nishimura in ACM Transactions
// on Modeling and Computer Simulation, vol. 8, no. 1,
// January 1998, pp. 3-30.
// And for the 64 bit version in the paper "Tables of 64-bit
// Mersenne Twisters" by Takuji Nishimura
// ************************************************************
// The original C code contained the following notice:
// Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
......@@ -18,69 +20,83 @@
implementation module MersenneTwister
import
StdArray, StdBool, StdEnum, StdDebug,
StdInt, StdList, StdMisc, StdReal, StdString
import StdArray,StdBool,StdEnum,StdDebug,StdInt,StdList,StdMisc,StdReal,StdString
n :== 624
m :== 397
N :== IF_INT_64_OR_32 312 624
M :== IF_INT_64_OR_32 156 397
a :== 0x9908b0df
uppermask :== 0x80000000
lowermask :== 0x7fffffff
Nm1 :== IF_INT_64_OR_32 311 623
NmM :== IF_INT_64_OR_32 156 227
MmN :== IF_INT_64_OR_32 -156 -227
mask_b :== 0x9d2c5680
mask_c :== 0xefc60000
A :== IF_INT_64_OR_32 0xB5026F5AA96619E9 0x9908b0df
uppermask :== IF_INT_64_OR_32 0xFFFFFFFF80000000 0x80000000
lowermask :== 0x7FFFFFFF
u :== 11
s :== 7
t :== 15
l :== 18
Mask_b :== IF_INT_64_OR_32 0x71D67FFFEDA60000 0x9d2c5680
Mask_c :== IF_INT_64_OR_32 0xFFF7EEE000000000 0xefc60000
shiftRight m n :== (m >> n) bitand (0x7FFFFFFF >> (n-1))
u :== IF_INT_64_OR_32 29 11
s :== IF_INT_64_OR_32 17 7
t :== IF_INT_64_OR_32 37 15
l :== IF_INT_64_OR_32 43 18
initrand :: Int *{# Int} -> *{# Int}
shiftRight :: !Int !Int -> Int;
shiftRight m n = code inline {
shiftrU
}
initrand :: Int *{#Int} -> *{#Int}
initrand seed mt
| seed==0 = abort "initrand: seed must not be 0\n"
= init_i 0 { mt & [0] = seed bitand 0xffffffff}
| seed==0
= abort "initrand: seed must not be 0\n"
= init_i 1 {mt & [0] = seed} seed
where
init_i i mt =: { [i] = mti}
| ii==n = mt
= init_i ii { mt & [ii] = (69069 * mti) bitand 0xffffffff}
where
ii = inc i
init_i :: !Int !*{#Int} !Int -> *{#Int}
init_i i mt n
| i<N
# n = IF_INT_64_OR_32 (n * 69069) ((n bitxor (shiftRight n 62)) * 6364136223846793005 + i);
= init_i (i+1) {mt & [i] = n} n
= mt
genRandInt :: Int -> [Int]
genRandInt seed
#! first_mt = initrand seed (createArray n 0)
#! first_mt = initrand seed (createArray N 0)
= gr_mti 0 (new_mti 0 first_mt)
where
gr_mti mti mt
| mti==n = gr_mti 0 (new_mti 0 mt)
#! (y, mt) = mt![mti]
#! y = y bitxor (shiftRight y u)
#! y = y bitxor ((y << s) bitand mask_b)
#! y = y bitxor ((y << t) bitand mask_c)
#! y = y bitxor (shiftRight y l)
= [ y : gr_mti (inc mti) mt]
new_mti :: Int *{# Int} -> *{# Int}
| mti==N
= gr_mti 0 (new_mti 0 mt)
#! y = mt.[mti]
# y = y bitxor (IF_INT_64_OR_32 ((shiftRight y u) bitand 0x5555555555555555) (shiftRight y u))
# y = y bitxor ((y << s) bitand Mask_b)
# y = y bitxor ((y << t) bitand Mask_c)
#! y = y bitxor (shiftRight y l)
= [y : gr_mti (inc mti) mt]
new_mti :: Int *{#Int} -> *{#Int}
new_mti k mt
| k==n = mt
#! (mtk, mt) = mt![k]
#! (mtkk, mt) = mt![(inc k) rem n]
#! y = (mtk bitand uppermask) bitor (mtkk bitand lowermask)
#! (mtind, mt) = mt![ind k]
#! mttmp = mtind bitxor (shiftRight y 1) bitxor ((y bitand 0x1) * a)
= new_mti (inc k) { mt & [k] = mttmp}
where
ind :: Int -> Int
ind k
| k<n-m = k+m
= k+m-n
| k<NmM
#! y = (mt.[k] bitand uppermask) bitor (mt.[k+1] bitand lowermask)
#! mttmp = mt.[k+M] bitxor (shiftRight y 1) bitxor ((~(y bitand 1)) bitand A)
= new_mti (inc k) {mt & [k] = mttmp}
= new_mti2 k mt
new_mti2 :: Int *{#Int} -> *{#Int}
new_mti2 k mt
| k<Nm1
#! y = (mt.[k] bitand uppermask) bitor (mt.[k+1] bitand lowermask)
#! mttmp = mt.[k+MmN] bitxor (shiftRight y 1) bitxor ((~(y bitand 1)) bitand A)
= new_mti2 (inc k) {mt & [k] = mttmp}
#! y = (mt.[k] bitand uppermask) bitor (mt.[0] bitand lowermask)
#! mttmp = mt.[k+MmN] bitxor (shiftRight y 1) bitxor ((~(y bitand 1)) bitand A)
= {mt & [k] = mttmp}
uint2Real i :== if (i >= 0) (toReal i) (4294967296.0 + toReal i)
genRandReal :: Int -> [Real]
genRandReal seed = map (\y -> (uint2Real y) / 4294967295.0) (genRandInt seed)
genRandReal seed
= [IF_INT_64_OR_32
((toReal (shiftRight y 11)) / 9007199254740991.0)
((uint2Real y) / 4294967295.0)
\\ y<-genRandInt seed]
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment