The Assembler version got through 1TB PractRand analysis with a small anamoly during the 256MB test. So, passed with flying colours.
Below is the source RndMT.bas but
here is a link, to save you the trouble, to RndMT.zip which includes a folder RndMT which, in turn includes, RndMT.bi and libRndMT.a.
I have limited experience with building libraries but so far the code seems to take a performance hit. Not so with RndMT - the library is as fast as the raw source code.
Here is a copy of the bi file:
Code: Select all
Declare Function RndMTS() As Single
Declare Function RndMTD() As Double
Declare Function RndMTR( As Long, As Long ) As Long
Declare Function RndMTUL() As ULong
Declare Sub RandomizeRndMT( As String )
#inclib "RndMT"
So, we have Single, Double, Range and a new one UL. There is probably not much use for RndMTUL but since the generator churns out 32 bit integers for a living it seemed a rightful addition.
RandomizeRndMT( "" ) ==> System seeding
RandomizeRndMT( "Some string" ) ==> User seeding
Example usage:
Code: Select all
#Include Once "RndMT.bi"
RandomizeRndMT( "" )
Print RndMTS
Print RndMTD
Print RndMTUL
Print RndMTR(1,100)
Print "Done"
Sleep
Code: Select all
RndMT FreeBasic
Entry points 2^19937 - 1 2^32
Speed 95.3MHz 85.8MHz (i7-3770K 3.50GHz)
Don't forget, Windows Vista or later. However, you have the source code so you can do whatever you like with it.
Have fun.
RndMT.bas
Code: Select all
' The C-program for MT19937, with initialization improved 2002/1/26
' was coded by Takuji Nishimura and Makoto Matsumoto.
'
' Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
' All rights reserved.
'
' Translated to FreeBASIC by Greg Lyon - April 2006
'
' UInteger conversion to Single and Double, Range and entry point
' to Mersenne Twister maximised by David Roberts - March 2017
Const N As uInteger = 624
Const M As uInteger = 397
Const MATRIX_A As uInteger = &H9908b0dfUL ' constant vector a
Const UPPER_MASK As uInteger = &H80000000UL ' most significant w-r bits
Const LOWER_MASK As uInteger = &H7fffffffUL ' least significant r bits
Dim Shared mt(N) As uInteger ' the array For the state vector
Dim Shared mti As Long = N + 1 ' mti == N+1 means mt[N] Is Not initialized
' initializes mt[N] with a seed
Private Sub init_genrand(ByVal s As uInteger)
mt(0)= s And &HffffffffUL
For mti = 1 To N-1
mt(mti) = ( 1812433253UL * ( mt(mti-1) Xor ( mt(mti-1) Shr 30 ) ) + mti )
' See Knuth TAOCP Vol2. 3rd Ed. P.106 For multiplier.
' In the previous versions, MSBs of the seed affect
' only MSBs of the array mt[].
' 2002/01/09 modified by Makoto Matsumoto
mt(mti) = mt(mti) And &HffffffffUL
' For >32 Bit machines
Next mti
End Sub
'-------------------------------------------------------
' Adaptation of code by Sterling Stuart Stein
Private Function genrand_int32() As UInteger
Static As UInteger Ptr ptrMT
Static mag01(1) As UInteger = {&h00000000, &h9908b0df}
Static As UInteger Ptr ptrmag01
Static As Long Counter = 0
ptrMT = @mt(0)
ptrmag01 = @mag01(0)
Asm
mov edi, dword ptr [ptrMT]
cmp dword ptr [mti], N
jb XB
mov esi, edi
SB:
mov eax,[esi]
and eax, UPPER_MASK
mov ebx,[esi+4]
and ebx, LOWER_MASK
or eax,ebx
mov ecx,eax
shr eax,1
mov edx,esi
add edx,(M*4)
mov ebx,[edx]
xor eax,ebx
and ecx,1
mov ebx, dword ptr [ptrmag01]
xor eax, [ebx + ecx*4]
mov [esi],eax
add esi,4
inc dword ptr [Counter]
cmp dword ptr [Counter],(N-M)
jnz SB
UB:
mov eax,[esi]
and eax,UPPER_MASK
mov ebx,[esi+4]
and ebx,LOWER_MASK
or eax,ebx
mov ecx,eax
shr eax,1
mov edx,esi
add edx,((M-N)*4)
mov ebx,[edx]
xor eax,ebx
and ecx,1
mov ebx, dword ptr [ptrmag01]
xor eax, [ebx + ecx*4]
mov [esi],eax
add esi,4
inc dword ptr [Counter]
cmp dword ptr [Counter],(N-1)
jnz UB
mov dword ptr [Counter], 0
' Added by David Roberts
'y = (mt((N - 1)) and UPPER_MASK) or (mt(0) and LOWER_MASK)
'---------------------------
mov eax, [esi]
and eax,UPPER_MASK
mov ebx, [edi]
and ebx, LOWER_MASK
or eax, ebx
mov ecx,eax
shr eax, 1
'---------------------------
mov edx,edi
add edx,(M-1)*4
mov ebx,[edx]
xor eax,ebx
and ecx,1
mov ebx, dword ptr [ptrmag01]
xor eax, [ebx + ecx*4]
mov [esi],eax
mov dword ptr [mti], 0
XB:
mov esi,edi
mov eax, dword ptr [mti]
inc dword ptr [mti]
shl eax,2
add esi,eax
mov eax,[esi]
mov ebx,eax
shr eax,11
xor ebx,eax
mov eax,ebx
shl eax,7
and eax,&H9d2c5680
xor ebx,eax
mov eax,ebx
shl eax,15
and eax,&Hefc60000
xor ebx,eax
mov eax,ebx
shr eax,18
xor eax,ebx
mov dword ptr [Function], eax
End Asm
End Function
'-------------------------------------------------------
' generates a Random number On [0,1)-real-interval
Public Function RndMTS() As Single
'Function genrand_real2() As Single
Static As UInteger TempVar
TempVar = genrand_int32()
' Adaptation of code by Wilbert @ PureBASIC forums
Asm
mov eax, dword ptr [TempVar]
movd xmm0, eax
psrlq xmm0, 9
mov eax, 1
cvtsi2ss xmm1, eax
por xmm0, xmm1
subss xmm0, xmm1
movd [Function], xmm0
End Asm
'Function = genrand_int32() * 1.0#/4294967296#
' divided by 2^32
End Function
'-------------------------------------------------------
' generates a random number on [0,1) with 53-bit resolution
'Function genrand_res53() As Double
Public Function RndMTD() As Double
Static As UInteger TempVar1, TempVar2
TempVar1 = genrand_int32()
TempVar2 = genrand_int32()
' Adaptation of code by Wilbert @ PureBASIC forums
Asm
mov eax, dword ptr [TempVar1]
movd xmm0, eax
mov eax, dword ptr [TempVar2]
movd xmm1, eax
punpckldq xmm0, xmm1
psrlq xmm0, 12
mov eax, 1
cvtsi2sd xmm1, eax
por xmm0, xmm1
subsd xmm0, xmm1
movq [Function], xmm0
End Asm
'Dim a As uInteger = genrand_int32() Shr 5
'Dim b As uInteger = genrand_int32() Shr 6
'Function = (a * 67108864.0# + b) * (1.0# / 9007199254740992.0#)
End Function
'-------------------------------------------------------
'The following code added by David Roberts
#include once "windows.bi"
#include once "win/bcrypt.bi"
#inclib "bcrypt"
Declare Function MyRandomInt Lib "Advapi32.dll" Alias "SystemFunction036" _
( RandomBuffer As Any Ptr, RandomBufferLength As ULong ) As Byte
Private Function SHA512( sText as String ) as String
Dim as Byte Ptr phalg, phhash
Dim as ULong lhashlength, lresult
Dim sbinhash As String
BCryptOpenAlgorithmProvider Varptr(phalg), wstr("SHA512"), 0, 0 ' we want phalg
BCryptCreateHash phalg, @phhash, null, 0, 0, 0, 0 ' we want phhash
BCryptHashData( phhash, Strptr( stext), len( stext ), 0 )
BCryptGetProperty phalg, bcrypt_hash_length, Cast ( Puchar, @lhashlength ), 4, @lresult, 0
sbinhash = String$( lhashlength, 0 )
BCryptFinishHash phhash, Strptr( sbinhash ), lhashlength, 0
BCryptDestroyHash phhash
BCryptCloseAlgorithmProvider phalg, 0
Return sbinhash
End Function
'-------------------------------------------------------
' Entry point to RndMT.bas. Next statement will be executed
' Initialize state vector in case RandomizeRndMT is not employed
init_genrand(5489UL)
'-------------------------------------------------------
Public Sub RandomizeRndMT( seedString As String )
Dim BinHash As String * 64
Dim As Ulong Ptr BinHashPtr, BaseBinHashPtr
Dim arrayPtr as Byte Ptr
Dim As Long i, j
If seedString <> "" Then ' User seed
BinHashPtr = Cptr(Ulong Ptr, StrPtr( BinHash ))
BaseBinHashPtr = BinHashPtr
BinHash = SHA512( seedString )
' j = 0 to 38 x i = 0 to 15 ==> 39 x 16 = 624 ( = N )
For j = 0 to 38
BinHash = SHA512( BinHash ) ' Hash the hash
BinHashPtr = BaseBinHashPtr
For i = 0 to 15
mt(j*16 + i) = PeeK( Ulong, BinHashPtr ) ' Populate state vector
BinhashPtr += 1
Next
Next
Else ' System seed
arrayPtr = Cast( Byte Ptr, @mt(0) )
MyRandomInt( arrayPtr, 4*N ) ' Populate state vector
End If
mti = N ' Signal state vector needs to be 'refurbished'
End Sub
'-------------------------------------------------------
Public Function RndMTR( One As Long, Two As Long ) As Long
Static As UInteger TempVar
TempVar = genrand_int32()
' Adaptation of code by John Gleason @ PowerBASIC forums
Asm
mov edx, dword ptr [TempVar]
mov ecx, dword ptr [One]
mov eax, dword ptr [Two]
cmp ecx, eax
jl Now1LT2
xchg eax, ecx
Now1LT2:
Sub eax, ecx
inc eax
jz doTheRnd
mul edx
Add edx, ecx
doTheRnd:
mov dword ptr [Function], edx
End Asm
End Function
'-------------------------------------------------------
Public Function RndMTUL() As ULong
Function = genrand_int32()
End Function
'-------------------------------------------------------
PractRand results
Code: Select all
Microsoft Windows [Version 10.0.14393]
(c) 2016 Microsoft Corporation. All rights reserved.
C:\Users\deltarho>F:
F:\>cd msvc12_32bit
F:\msvc12_32bit>My_Rng | RNG_test stdin32
RNG = RNG_stdin32, PractRand version 0.91, seed = 0xd5bb4d5f
test set = normal, folding = standard (32 bit)
rng=RNG_stdin32, seed=0xd5bb4d5f
length= 128 megabytes (2^27 bytes), time= 3.0 seconds
no anomalies in 117 test result(s)
rng=RNG_stdin32, seed=0xd5bb4d5f
length= 256 megabytes (2^28 bytes), time= 6.3 seconds
Test Name Raw Processed Evaluation
[Low1/32]DC6-9x1Bytes-1 R= +6.3 p = 3.1e-3 unusual
...and 123 test result(s) without anomalies
rng=RNG_stdin32, seed=0xd5bb4d5f
length= 512 megabytes (2^29 bytes), time= 12.5 seconds
no anomalies in 132 test result(s)
rng=RNG_stdin32, seed=0xd5bb4d5f
length= 1 gigabyte (2^30 bytes), time= 24.8 seconds
no anomalies in 141 test result(s)
rng=RNG_stdin32, seed=0xd5bb4d5f
length= 2 gigabytes (2^31 bytes), time= 48.7 seconds
no anomalies in 148 test result(s)
rng=RNG_stdin32, seed=0xd5bb4d5f
length= 4 gigabytes (2^32 bytes), time= 95.6 seconds
no anomalies in 156 test result(s)
rng=RNG_stdin32, seed=0xd5bb4d5f
length= 8 gigabytes (2^33 bytes), time= 190 seconds
no anomalies in 165 test result(s)
rng=RNG_stdin32, seed=0xd5bb4d5f
length= 16 gigabytes (2^34 bytes), time= 379 seconds
no anomalies in 172 test result(s)
rng=RNG_stdin32, seed=0xd5bb4d5f
length= 32 gigabytes (2^35 bytes), time= 751 seconds
no anomalies in 180 test result(s)
rng=RNG_stdin32, seed=0xd5bb4d5f
length= 64 gigabytes (2^36 bytes), time= 1509 seconds
no anomalies in 189 test result(s)
rng=RNG_stdin32, seed=0xd5bb4d5f
length= 128 gigabytes (2^37 bytes), time= 3015 seconds
no anomalies in 196 test result(s)
rng=RNG_stdin32, seed=0xd5bb4d5f
length= 256 gigabytes (2^38 bytes), time= 6001 seconds
no anomalies in 204 test result(s)
rng=RNG_stdin32, seed=0xd5bb4d5f
length= 512 gigabytes (2^39 bytes), time= 12037 seconds
no anomalies in 213 test result(s)
rng=RNG_stdin32, seed=0xd5bb4d5f
length= 1 terabyte (2^40 bytes), time= 24056 seconds
no anomalies in 220 test result(s)