Mersenne Twister (erroneous, don't use)

[Hades]

Hi!

As I don't have enough information about PB's PRNG, I have converted the Mersenne Twister to PB. The generated numbers differ from the original, though. Hopefully just because PB doesn't have unsigned longs, but I can't be sure.
I have done some tests, and everything looks fine. But maybe some of you have different tests, and can spot errors?

Here, PBs Random() is 3-4 times faster than the MT-version, but the Mersenne Twister has a very long period (2^19937-1, time until it repeats), and generates 'better' random numbers than most other PRNGs (I have not enough information to compare it to PBs, though).

And with ASM, it's speed could be improoved by a large margin, so there probably wouldn't be to much difference to PBs PRNG.

Code: Select all

;             Mersenne Twister 
;    pseudo-random number generator (PRNG)
;
; PureBasic conversion by Hades, May 2006, 
; of... 
;
; /*
; A C-program for MT19937, with initialization improved 2002/1/26.
; Coded by Takuji Nishimura and Makoto Matsumoto.
; 
; Before using, initialize the state by using init_genrand(seed)  
; or init_by_array(init_key, key_length).
; 
; Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
; All rights reserved.                           
; Copyright (C) 2005, Mutsuo Saito,
; All rights reserved.                         
; 
; Redistribution and use in source and binary forms, with or without
; modification, are permitted provided that the following conditions
; are met:
; 
;   1.  Redistributions of source code must retain the above copyright
;       notice, this list of conditions and the following disclaimer.
; 
;   2.  Redistributions in binary form must reproduce the above copyright
;       notice, this list of conditions and the following disclaimer in the
;       documentation and/or other materials provided with the distribution.
; 
;   3.  The names of its contributors may not be used to endorse or promote 
;       products derived from this software without specific prior written 
;       permission.
; 
; THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
; "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
; LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
; A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
; CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
; EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
; PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
; LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
; NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
; 
; 
; Any feedback is very welcome.
; http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
; email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
; */
     



; Period parameters 
#MT_N = 624
#MT_M = 397
#MT_MatrixA = $9908B0DF     ; constant vector a
#MT_UpperMask = $80000000   ; most significant w-r bits
#MT_LowerMask = $7FFFFFFF   ; least significant r bits

Global Dim MT_mt.l(#MT_N - 1)   ; the array for the state vector
Global MT_mti.l = #MT_N + 1     ; mti==N+1 means mt[N] is not initialized

Global Dim MT_mag01(1)  ; MT_mag01[x] = x * #MT_MatrixA  for x=0,1
MT_mag01(0) = 0
MT_mag01(1) = #MT_MatrixA



Procedure MT_init_genrand(Seed.l) ; initializes the Mersenne Twister with a seed (-2147483648 <= Seed <= 2147483647)
  MT_mt(0) = Seed
  For MT_mti = 1 To #MT_N-1
    MT_mt(MT_mti) = 1812433253 * (MT_mt(MT_mti-1) ! (((MT_mt(MT_mti-1) >> 1) & #MT_LowerMask) >> 29)) + MT_mti
  Next
EndProcedure


; initialize by an array with array-length
; init_key is the array for initializing keys
; key_length is its length 
Procedure MT_init_by_array(*init_key.l, key_length.l)
  Protected i.l, j.l, k.l, n.l
  MT_init_genrand(19650218)
  i = 1 : j = 0
  k = key_length
  If k > #MT_N 
    k = #MT_N
  EndIf
  For n = k To 0 Step -1
    MT_mt(i) = (MT_mt(i) ! ((((MT_mt(i-1) >> 1) & #MT_LowerMask) >> 29) * 1664525)) + PeekL(*init_key + 4 * j) + j
    i + 1 : j + 1
    If i >= #MT_N 
      MT_mt(0) = MT_mt(#MT_N-1)
      i = 1
    EndIf
    If j >= key_length
      j = 0
    EndIf
  Next
  For k = #MT_N - 1 To 0 Step -1
    MT_mt(i) = (MT_mt(i) ! (MT_mt(i-1) ! (((MT_mt(i-1) >> 1) & #MT_LowerMask) >> 29) * 1566083941)) - i
    i + 1
    If i >= #MT_N
      MT_mt(0) = MT_mt(#MT_N - 1)
      i = 1
    EndIf
  Next
  MT_mt(0) = #MT_UpperMask
EndProcedure


Procedure.l MT_genrand_int32() ; Returns a random long (0 .. $7FFFFFFF)
  Protected y.l, kk.l, yh.l
  If MT_mti >= #MT_N
    If MT_mti = #MT_N + 1   ; if MT_init_genrand() has not been called,
      MT_init_genrand(5489) ; a default initial seed is used
    EndIf
    For kk = 0 To (#MT_N - #MT_M) - 1
      y = (MT_mt(kk) & #MT_UpperMask) | (MT_mt(kk+1) & #MT_LowerMask)
      MT_mt(kk) = MT_mt(kk+#MT_M) ! ((y >> 1) & #MT_LowerMask) ! MT_mag01(y & 1)
    Next
    For kk = #MT_N - #MT_M To #MT_N - 2
      y = (MT_mt(kk) & #MT_UpperMask) | (MT_mt(kk+1) & #MT_LowerMask)
      MT_mt(kk) = MT_mt(kk+(#MT_M - #MT_N)) ! ((y >> 1) & #MT_LowerMask) ! MT_mag01(y & 1)
    Next
    y = (MT_mt(#MT_N-1) & #MT_UpperMask) | (MT_mt(0) & #MT_LowerMask)
    MT_mt(#MT_N-1) = MT_mt(#MT_M-1) ! ((y >> 1) & #MT_LowerMask) ! MT_mag01(y & 1)
    MT_mti = 0
  EndIf
  y = MT_mt(MT_mti)
  MT_mti + 1
  ; Tempering
  y = y ! (((y >> 1) & #MT_LowerMask) >> 10)
  yh = y << 6
  If yh And $40000000
    y = y ! (((yh << 1) | #MT_UpperMask) & $9D2C5680)
  Else
    y = y ! ((yh << 1) & $9D2C5680)
  EndIf
  yh = y << 14
  If yh And $40000000
    y = y ! (((yh << 1) | #MT_UpperMask) & $EFC60000)
  Else
    y = y ! ((yh << 1) & $EFC60000)
  EndIf
  y = y ! (((y >> 1) & #MT_LowerMask) >> 17)
  
  ProcedureReturn y & #MT_LowerMask
EndProcedure


Procedure.f MT_Rnd() ; Returns a random float (0.0 .. 1.0)
  ProcedureReturn MT_genrand_int32() * (1.0 / 2147483647.0)
EndProcedure

Procedure.l MT_Random(Max.l) ; Returns a random long (0 .. Max, 1 <= Max <= 2147483647) 
  If Max <= 0
    ProcedureReturn 0
  EndIf
  ProcedureReturn MT_genrand_int32() % (Max + 1)
EndProcedure

[Pupil]

It's probably the '>>' operator that make you get different numbers in PB. (PB uses arithmetic shift right for this operator, which means that the sign bit is preserved i.e. one's are shifted in from the left)

[yoxola]

Nice idea to use Mersenne Twister

I use milisecs as seed to generate random number,
and that's not random enough...

[Dare2]

Nice.

Are duplicates within a largish range (on same seed) of concern to you?

If so, with this test:

Code: Select all

Dim res.l(52000)
MT_init_genrand($7FFFFFFF)
For i=1 To 52000
  res(i) = MT_Random($7FFFFFFF)
Next
SortArray(res(),0)
last.l = $7FFFFFFF
ctr=0
For i = 1 To 52000
  If res(i)=last
  Debug Str(i)+":"+Str(last)+"="+Str(res(i))
    ctr+1
  EndIf
  last=res(i)
Next
Debug ctr

I get a duplicate. With PureBasic, I get 2 using same seed and random call.

Seeded with $FFFFFFFFF get none with yours, and Pure gives 1.

[Hades]

@Pupil

Because of that I have changed all :
'y >> x'
to :
'((y >> 1) & $7FFFFFFF) >> (x-1)'
and all :
'y << x'
to :
' yh = y << (x-1)
If yh And $40000000
y = ((yh << 1) | $80000000)
Else
y = (yh << 1)
EndIf
'
Thats the reason for it being slow and ugly in PB. :roll:
But inside the init... procedures there is multiplication, and that makes a difference if unsigned. But it shoudn't break the MT.

@Dare2

I have to create probably thousands of numbers for a single seed, so if the first one or two are duplicates it will be no problem. However, if most or all numbers for different seeds would be the same, it would be a big problem for me. With the MT this shouldn't happen, if I have done it right.

Btw, thanks for testing.

[Hades]

Hi,

mskuma pointed me to ComScire RNGmeter (lower half of that site). A nice program to test random number generators.

It takes a file with random numbers and runs some tests on them. The more tests are passed, the higher the rating.

My implementation of the MT didn't fail for 100 million numbers and got a rating of 29,6+ ('+' means higher than). I haven't tested any further yet.

In comparsion using PBs Random($FFFFFFFF) fails for less than 10000 numbers and gets a rating of 14.0- ('-' means lower than).


So this MT implementation should be save to use for most cases except cryptography and online gambling.


Have fun,

Hades

[mskuma]

In comparsion using PBs Random($FFFFFFFF) fails for less than 10000 numbers and gets a rating of 14.0- ('-' means lower than).

I'd say 10,000 isn't enough for a decent test, and probably the test should be done using the same quantity of numbers. I tried 100 million numbers using PB's RNG and got a score of 27.6+. Also I got a 14 (ambivalent, no plus or minus) when testing for 10,000 MT numbers and PB numbers.

Cheers.

[Hades]

Strange...

The test runs until it fails. So if it fails (detects errors or regularities) after less than 10,000 numbers it doesn't matter if you try millions of numbers. It will still fail after less than 10,000, and will get the same score.

I did a lot of tests with different numers, and it was always the same: the MT didn't fail, PBs PRNG always failed.
I just did some more tests, and PBs PRNG didn't failed for 2500 numbers and got a score of 0.0+. 3000 numbers and everything above fails, and gives a score of 14.0- .


Edit: Ok, it was my mistake.
PBs PRNG doesn't create 32 random bits (probably 31), and the test program recognizes that. When I'm creating the numbers using '(Random($FFFF) << 16) + Random($FFFF)' PBs PRNG doesn't fail for 100 Million numbers, too. So PBs PRNG is all fine.

Thanks for rechecking, mskuma.

Edit: mskuma did some further tests, and showed that my port of the Mersenne Twister doesn't behave like expected. I don't know how serious the problem is, and how to solve it, so better don't use it.
Also he sent me a link to a dll version of the MT, so it doesn't make much sense for me to continue work on 'my own' version.

Bye,

Hades