[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Octave-bug-tracker] [bug #42557] rand seed for twister algorithm
From: |
Markus Bergholz |
Subject: |
[Octave-bug-tracker] [bug #42557] rand seed for twister algorithm |
Date: |
Sun, 15 Jun 2014 08:28:32 +0000 |
User-agent: |
Mozilla/5.0 (X11; Linux x86_64; rv:30.0) Gecko/20100101 Firefox/30.0 |
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/
- [Octave-bug-tracker] [bug #42557] rand seed for twister algorithm,
Markus Bergholz <=