Re: Same sequence of random number generators than under Matlab?

[Francesco Potortì]

>> x_{k+1} = ax_k (mod m) 
>
>
>> with a,m = 16807, 2^31-1 
>
>I tried to implement the algorithm (see attached file) but it doesn't
>seem to work (I get sequences containing 0). Can somebody help?

octave> format long
octave> x=1
x =  1
octave> x=mod(16807*x,2^31-1)
x =  16807
octave> x=mod(16807*x,2^31-1)
x =  282475249
octave> x=mod(16807*x,2^31-1)
x =  1622650073
octave> x=mod(16807*x,2^31-1)
x =  984943658
octave> x=mod(16807*x,2^31-1)
x =  1144108930
octave> x=mod(16807*x,2^31-1)
x =  470211272
octave> x=mod(16807*x,2^31-1)
x =  101027544
octave> x=mod(16807*x,2^31-1)
x =  1457850878
octave> x=mod(16807*x,2^31-1)
x =  1458777923
octave> tic; x(1)=1; for i=2:1000000; x(i)=mod(16807*x(i-1),2^31-1); endfor; toc
Elapsed time is 13.190524101 seconds.
octave> statistics(x)
ans =

 Columns 1 through 3:

   1.00000000000000e+00   5.37129909750000e+08   1.07464864950000e+09

 Columns 4 through 6:

   1.61026491475000e+09   2.14748353100000e+09   1.07380514916780e+09

 Columns 7 through 9:

   6.19606696291176e+08  -5.71388955027976e-04  -1.19893487840676e+00

-- 
Francesco Potortì (ricercatore)        Voice:  +39.050.315.3058 (op.2111)
ISTI - Area della ricerca CNR          Mobile: +39.348.8283.107
via G. Moruzzi 1, I-56124 Pisa         Fax:    +39.050.315.2040  
(entrance 20, 1st floor, room C71)     Web:    http://fly.isti.cnr.it