Re: redefining rand and srand

Hi,

"Phil Tomson" <[e-mail]> wrote in message
news:[e-mail]...
> In article <[e-mail]>,
> >
> >A couple of questions:
> >1) 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.
> >2) 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] = 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 (a*67108864.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