# bug #42557: rand seed for twister algorithm differs from Matlab

 Submitted by: Markus Bergholz Submitted on: Sun 15 Jun 2014 08:28:31 AM UTC Category: Octave Function Severity: 3 - Normal Priority: 5 - Normal Item Group: Matlab Compatibility Status: Confirmed Assigned to: None Originator Name: Open/Closed:  * Open Release:  * dev Operating System:  * Any
* Mandatory Fields

Comment Type & Canned Response:

## Discussion

Thu 19 Nov 2015 05:34:58 PM UTC, comment #7:

Actually, scrap that last comment, 0 is seed 5489 in Matlab's latest RNG: http://www.mathworks.com/matlabcentral/answers/28216-how-to-duplicate-rand-twister-0-in-current-matlab#answer_36524

So only the normal distribution discrepancy remains.

Anonymous
Thu 19 Nov 2015 03:48:12 PM UTC, comment #6:

The flipping of the last element of the state vector reference was supposed to link to bug #44046

Anonymous
Thu 19 Nov 2015 02:04:17 PM UTC, comment #5:

This doesn't seem to hold for randn, tried with

Octave 4.0.0:

>> randn('twister', twister_seed(5489));
>> randn(1,5)

ans =

1.362259   0.071959   0.203598  -0.096973   1.324298

Matlab R2015a:

>> rng(5489, 'twister');
>> randn(1,5)

ans =

0.5377    1.8339   -2.2588    0.8622    0.3188

Also tried flipping the last state vector value as decribed (in both Matlab and Octave), randn still produces different streams.

Matlab also seems to use the Ziggurat method, so can't figure out where the difference is.

>> RandStream.getGlobalStream

ans =

mt19937ar random stream (current global stream)
Seed: 5489
NormalTransform: Ziggurat

Oddly enough, Matlab always seems to produce the same randn sequence no matter what the seed number is, i.e. rng(0, 'twister') and rng(5489, 'twister') generate the same 625 state vector when queried with s = rng

Anonymous
Sat 28 Jun 2014 05:15:20 AM UTC, comment #4:

>> rand ('twister', [1 2 3]);
Error using rand
State must be a double scalar or the output of RAND('twister').

>> rng([1 2 3])
Error using rng (line 128)
First input must be a nonnegative integer seed less than 2^32, 'shuffle', 'default', or generator settings captured previously using S = RNG.

Markus Bergholz <markuman>
Thu 26 Jun 2014 03:41:21 PM UTC, comment #3:

What does Matlab do if the seed is a vector, rather than a scalar?

rand ('twister', [1 2 3]);
rand

Rik <rik5>
Fri 20 Jun 2014 05:37:20 PM UTC, comment #2:

and finaly

octave:8> rand('twister',twister_seed(5489))
octave:9> rand(5,1)
ans =

0.814723686393179
0.905791937075619
0.126986816293506
0.913375856139019
0.63235924622541

seems to be the default seed from matlab...

Markus Bergholz <markuman>
Mon 16 Jun 2014 12:25:10 AM UTC, comment #1:

I've marked the bug as confirmed and changed the category to "Matlab Compatibility".

Rik <rik5>
Sun 15 Jun 2014 08:28:31 AM UTC, original submission:

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

Markus Bergholz <markuman>

## Attached Files

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Comment:

Attached Files
file #31549:  twister_seed.m added by markuman (456B - text/x-objcsrc)

## Dependencies

Depends on the following items: None found

Items that depend on this one

Digest:
bug dependencies.

Carbon-Copy List
• -email is unavailable- added by rik5 (Posted a comment)
• -email is unavailable- added by markuman (Submitted the item)
•

There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

Only project members can vote.