Hi,
“Phil Tomson” ptkwt@shell1.aracnet.com wrote in message
news:b0ivi80pgg@enews3.newsguy.com…
In article b0idhk0gj3@enews4.newsguy.com,
A couple of questions:
- how good is Ruby’s rand function (what’s the quality of randomness)?
I
suspec that it is just the same as C’s rand, so there won’t be a
difference.
- I’d like to plug in my own random number generator that would probably
produce higher-quality random numbers than the current rand does. I’m
pretty sure it should be doable and easy, like:
module Kernel
def rand(max=0)
#do random number generation magic here
end
end
Here’s what I did. At Robert Feldt’s suggestion I used his RandomR
package which uses the Mersenne Twister to generate random numbers (a much
better method than whatever is built-in, I’m sure). I then redefined rand
and srand like so:
require ‘random/mersenne_twister’
module Kernel
def rand(max=0)
if not defined? @__random
@__rng=Random::MersenneTwister.new Time.now.to_i
end
@__rng.rand(max)
end
def srand(seed)
@__rng=Random::MersenneTwister.new seed
end
end
How about use native Ruby code instead of extension module?
The following is my rough translation from C code(mt19937ar.c) at
http://www.math.keio.ac.jp/~matumoto/emt.html
···
=================================================================
Period parameters
N = 624
M = 397
MATRIX_A = 0x9908b0df # constant vector a
UPPER_MASK = 0x80000000 # most significant w-r bits
LOWER_MASK = 0x7fffffff # least significant r bits
static unsigned long mt[N]; # the array for the state vector
static int mti=N+1; # mti==N+1 means mt[N] is not initialized
$mt = Array.new(N)
$mti = N+1
initializes mt[N] with a seed
def init_genrand(s)
$mt[0]= s & 0xffffffff
for $mti in 1…N
$mt[$mti] =
(1812433253 * ($mt[$mti-1] ^ ($mt[$mti-1] >> 30)) + $mti)
# See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
# In the previous versions, MSBs of the seed affect
# only MSBs of the array mt.
# 2002/01/09 modified by Makoto Matsumoto
$mt[$mti] &= 0xffffffff
# for >32 bit machines
end
$mti = N
end
initialize by an array with array-length
init_key is the array for initializing keys
key_length is its length
def init_by_array(init_key, key_length)
init_genrand(19650218)
i=1
j=0
k = (N>key_length ? N : key_length)
k.downto(1) {
$mt[i] = ($mt[i] ^ (($mt[i-1] ^ ($mt[i-1] >> 30)) * 1664525)) +
init_key[j] + j # non linear
$mt[i] &= 0xffffffff # for WORDSIZE > 32 machines
i+=1
j+=1
if (i>=N)
$mt[0] = $mt[N-1]
i=1
end
if (j>=key_length)
j=0
end
}
(N-1).downto(1) {
$mt[i] = ($mt[i] ^ (($mt[i-1] ^ ($mt[i-1] >> 30)) * 1566083941)) - i
non linear
$mt[i] &= 0xffffffff # for WORDSIZE > 32 machines
i+=1
if (i>=N)
$mt[0] = $mt[N-1]
i=1
end
}
$mt[0] = 0x80000000 # MSB is 1; assuring non-zero initial array
end
generates a random number on [0,0xffffffff]-interval
def genrand_int32()
mag01=[0x0, MATRIX_A]
# mag01 = x * MATRIX_A for x=0,1
if ($mti >= N) # generate N words at one time
if ($mti == N+1) # if if init_genrand() has not been called,
init_genrand(5489) # a default initial seed is used
end
for kk in 0...N-M
y = ($mt[kk]&UPPER_MASK)|($mt[kk+1]&LOWER_MASK)
$mt[kk] = $mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1]
end
for kk in (N-M) ... (N-1)
y = ($mt[kk]&UPPER_MASK)|($mt[kk+1]&LOWER_MASK)
$mt[kk] = $mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1]
end
y = ($mt[N-1]&UPPER_MASK)|($mt[0]&LOWER_MASK)
$mt[N-1] = $mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1]
$mti = 0
end
y = $mt[$mti]
$mti+=1
y ^= (y >> 11)
y ^= (y << 7) & 0x9d2c5680
y ^= (y << 15) & 0xefc60000
y ^= (y >> 18)
return y
end
generates a random number on [0,0x7fffffff]-interval
def genrand_int31()
return (genrand_int32()>>1)
end
generates a random number on [0,1]-real-interval
def genrand_real1(void)
return genrand_int32()*(1.0/4294967295.0)
# divided by 2^32-1
end
generates a random number on [0,1)-real-interval
def genrand_real2()
return genrand_int32()*(1.0/4294967296.0)
# divided by 2^32
end
generates a random number on (0,1)-real-interval
def genrand_real3()
return (genrand_int32() + 0.5)*(1.0/4294967296.0)
# divided by 2^32
end
generates a random number on [0,1) with 53-bit resolution
def genrand_res53()
a=genrand_int32()>>5
b=genrand_int32()>>6
return (a67108864.0+b)(1.0/9007199254740992.0)
end
init=[0x123, 0x234, 0x345, 0x456]
length=4
init_by_array(init, length)
printf("1000 outputs of genrand_int32()\n")
for i in 0...1000
printf("%10u ", genrand_int32())
if (i%5==4)
printf("\n")
end
end
printf("\n1000 outputs of genrand_real2()\n")
for i in 0 ...1000
printf("%10.8f ", genrand_real2())
if (i%5==4)
printf("\n")
end
end
=====================================================
Park Heesob