[Octave-bug-tracker] [bug #42557] rand seed for twister algorithm

[Markus Bergholz]

URL:
  <http://savannah.gnu.org/bugs/?42557>

                 Summary: rand seed for twister algorithm
                 Project: GNU Octave
            Submitted by: markuman
            Submitted on: Sun 15 Jun 2014 08:28:31 AM GMT
                Category: Octave Function
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Incorrect Result
                  Status: None
             Assigned to: None
         Originator Name: 
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any
                 Release: 3.8.1
        Operating System: Any

    _______________________________________________________

Details:

Initial situation:  

Ruby, Python and Numpy using twister algorithm as default.
We can force this algorithm in Octave and Matlab too.
So, if you seed each software with the value 0, it should return the same
random numbers.


Octave 3.8.1

octave:1> rand('twister',0)
octave:2> rand
ans =    0.844421851525048



Python 3.4.1

>>> import random
>>> random.seed(0)
>>> random.random()
0.8444218515250481



Python 3.4.1 Numpy 1.8.1

>>> import numpy
>>> numpy.random.seed(0)
>>> numpy.random.random()
0.5488135039273248




Ruby 2.1.2

irb(main):001:0> prng = Random.new(0)
=> #<Random:0x00000001d469e8>
irb(main):002:0> prng.rand()
=> 0.5488135039273248



Matlab 2013a

>> rand('twister',0)
>> rand

ans =

    0.548813503927325



So Octave is handling the seed input wrong. Afaiu the seed is processed in
http://hg.octave.org/octave/file/03c2671493f9/liboctave/cruft/ranlib/getsd.f
due use_old_generators. But maybe I'm wrong. However, This behavior can be
fixed with another function.


function ret = twister_seed(SEED=0)

    ret = uint32(zeros(625,1));
    ret(1) = SEED;
    for N = 1:623
        ## initialize_generator
        # bit-xor (right shift by 30 bits)
        uint64(1812433253)*uint64(bitxor(ret(N),bitshift(ret(N),-30)))+N; #
has to be uint64, otherwise in 4th iteration hit maximum of uint32!
        ret(N+1) = uint32(bitand(ans,uint64(intmax('uint32')))); # untempered
numbers
    endfor
    ret(end) = 1;   

endfunction


This function needs to be translated into c++ (that's past my c++ skills!) and
need to be called in line 276 of libinterp/corefcn/rand.cc and returned to
ColumnVector s (afaiu).
http://hg.octave.org/octave/file/03c2671493f9/libinterp/corefcn/rand.cc#l276



Octave 3.8.1

octave:1> rand('twister',twister_seed) # notice: default seed is 0 in the
function
octave:2> rand
ans =    0.548813503927325
octave:3> rand
ans =    0.715189366372419
octave:4> rand
ans =    0.602763376071644
octave:5> rand('twister',twister_seed(42))
octave:6> rand
ans =    0.374540118847363
octave:7> rand
ans =    0.950714306409916
octave:8> rand
ans =    0.731993941811405




Matlab 2013a

>> rand('twister',0)
>> rand

ans =

   0.548813503927325

>> rand

ans =

   0.715189366372419

>> rand

ans =

   0.602763376071644

>> rand('twister',42)
>> rand

ans =

         0.374540118847362

>> rand

ans =

         0.950714306409916

>> rand

ans =

         0.731993941811405





Python 3.4.1 Numpy 1.8.1

>>> import numpy
>>> numpy.random.seed(0)
>>> numpy.random.random()
0.5488135039273248
>>> numpy.random.random()
0.7151893663724195
>>> numpy.random.random()
0.6027633760716439
>>> numpy.random.seed(42)
>>> numpy.random.random()
0.3745401188473625
>>> numpy.random.random()
0.9507143064099162
>>> numpy.random.random()
0.7319939418114051





    _______________________________________________________

File Attachments:


-------------------------------------------------------
Date: Sun 15 Jun 2014 08:28:31 AM GMT  Name: twister_seed.m  Size: 456B   By:
markuman

<http://savannah.gnu.org/bugs/download.php?file_id=31549>

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?42557>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/