Tweet
The Multiply-with-Carry (MWC) pseudo random number generator (PRNG) by George Marsaglia has a couple of issues: the first is the use of 2^32 as a base which is resolved when 2^32-1 is used; the second is with the method itself which is resolved by using what is called Complementary Multiply-with-Carry (CMWC).
The period of the original idea is just less than 2^64 which is the square of PB's RND at 2^32. On my machine doing nothing else but compute RND the cycle starts to repeat after only a couple of minutes and, of course, there are only 2^32 possible entry points into the cycle.
It is worth noting that when using Randomize Timer we are restricting the number of entry points to seconds in a day x 64, the usual resolution of the system clock, giving then only 55259600 or 2^22.4 entry points.
The period of both MWC and CMWC can be enhanced by using a lag method where we use initial values x(0), x(1), ..., x(r-1) and an initial value for the carry corresponding with x(r-1). The original idea then is a lag~1.
I have been posting slls of late but two function 'shells', CMWC( Opt OptionalLag As Long ) and CWMCSeed( Opt Text As String ), are presented here, so you can kick them around to suit your requirements.
The code presented here uses a CMWC base 2^32-1 with lag~r where r = 2, 4, 8, ..., 4096. With r = 4096 we have a period of 2^131086 compared with 2^19937 for the Mersenne Twister PRNG.
The code assumes at least Windows XP SP3 and the latest compilers but can be re-written to accomodate all versions of Windows and compilers but with more effort.
In Marsaglia's post here we have
Above that section of code is a MWC base 2^32 lag~256 which uses a multiplier 'a' of 809430660. The carrys generated will never exceed the multiplier and it is therefore preferable to use an intial carry less than the multipler. However, the choice in the above is the same as with the MWC example and the used value of c is 362436; which is greater than the multiplier 18782. The term c=(t>>32) is what we would use for a base 2^32. I wonder then if the CMWC was a copy and paste and not fully edited. I shall come back to c=(t>>32).
The principle of CMWC is x = b - 1 - x where b is the base and x, on the right hand side, is the number generated by the MWC. With b = 2^32-1 we get x = 2^32 - 2 - x. 2^32 - 2 is 0xfffffffe; the value used for r in 'return(Q[i]=r-x'. The x's generated by MWC base 2^32-1 will be in the range [0,2^32-2]. If we use x = 0 and x = 2^32-2 in the CMWC we get [2^32-2,0]. So, the above is using a base of 2^32-1.
The term x=t+c had me fooled for a while, I was expecting a modulo.
In general, j = k*(j\k) + j MOD k ie quotient + remainder
Rearranging, j MOD k = j - k*(j\k)
With k = 2^32-1 we get
j MOD (2^32-1) = j - (2^32-1)*(j\(2^32-1))
= j - 2^32*(j\(2^32-1)) + j\(2^32-1)
Substituting Marsaglia's variables we get
x = t MOD (2^32-1) = t - 2^32*(t\(2^32-1)) + t\(2^32-1)
= t - 2^32*c + c ; with c = t\(2^32-1)
= t + c
Since 2^32*c = 0 because c is 32 bit.
I had not seen that before - it is neat isn't it?
This further confirms a base of 2^32-1.
In the following code x = t Mod &H0FFFFFFFF is used rather than x = t + c because it was found to be marginally faster; we are using Register t As Ext. I should imagine that x = t + c will come into its own with an assembler version.
So, why is c=(t>>32) being used instead of c = t\&H0FFFFFFFF? ie t\(2^32-1).
Answers on a post card to ......
In the following code c = t\&H0FFFFFFFF is used.
I should add that every implementation seen, from using C/C++, PureBASIC and Java, port Marsaglia's post as written - no questions asked.
At the moment I don't understand 'if(x<c){x++;c++;}' but have accepted it. More thought is required. Perhaps someone can put me wise on that?
For byte stream testing purposes we need to map [0,2^32-2] into [0,2^32-1] otherwise the fourth byte of a dword dump will be in the range [0,254] and not [0,255].
In addition to George Marsaglia's DIEHARD tests I use ent.exe from Fourmilab.
Here are the results from ten tests of 15MB dumps; five done one day and another five done another day.
For a full description of the tests see Fourmilab's home page.
With the chi square, in particular, non-random behaviour will put us into the regions < 10% and > 90%. We only have two in the above and they would qualify as almost suspect. The more tests we do the more likely one would qualify as suspect. There is nothing unusual with that.
The conclusion here is that there is no evidence to contradict the null hypothesis that CMWC produces good quality pseudo random numbers.
Compare with these tests of 15MB dumps.
These were 15MB dumps from here; a quantum random number generator service - the reason why I am using 15MB dumps since the service gives dumps of 1MB, 15MB and 100MB.
That chi square of 97.13% puts us into the suspect region but it may be a while before we see that again.
Obviously, there are many ways we can populate Q() with the initial random seeds. In the following code the CryptGenRandom API is used. It is easy to implement and, by definition, does not require seeding itself.
Best practice would call the function CMWC once in an application's initial stages so as to initialize the generator. It has an optional lag parameter which defaults to lag~256 if not given or if an invalid lag is given.
On my machine I am getting about 10 million dwords a second, irrespective of which lag is used - a bit of assembler would help.
So, which lag to use? I had reservations about a small multipler but they are not borne out in tests. There are some scenarios where exceptionally large periods are required where even Mersenne Twister is left wanting. A lag of 1024 and above exceed Mersenne Twister. Devices limited in memory would benefit from low memory Q() usage. A lag~64 needs only 256 bytes for Q() and has a period of 2^2077. You decide.
Post any comments here. I have one of my own.
Example code:
***********************************************************************************
There are situations where we need to have the same sequence per instance of generation along the lines of 'Randomize [number]'. The function CMWC will not do that as it stands. The CryptGenRandom aspect will need to be replaced.
One way is presented with the function CMWCSeed along the lines of 'Randomize [alphanumeric]' which, of course, we do not have with PB.
The principle is to take an alphanumeric seed and hash it with, say, SHA256 to populate a lag~8 Q(). This has a period of 2^285. To put that into perspective 2^285/2^32 = 2^253 ~ 10^76 times the period of RND. A SHA512 approach will give ~ 10^162 times the period of RND.
The alphanumeric parameter is optional. If not given "The time has come the walrus said...." is used; and why not?
Obviously, there are many ways to populate Q() with a seeded approach - using a hash is an easy way.
Example code:
The period of the original idea is just less than 2^64 which is the square of PB's RND at 2^32. On my machine doing nothing else but compute RND the cycle starts to repeat after only a couple of minutes and, of course, there are only 2^32 possible entry points into the cycle.
It is worth noting that when using Randomize Timer we are restricting the number of entry points to seconds in a day x 64, the usual resolution of the system clock, giving then only 55259600 or 2^22.4 entry points.
The period of both MWC and CMWC can be enhanced by using a lag method where we use initial values x(0), x(1), ..., x(r-1) and an initial value for the carry corresponding with x(r-1). The original idea then is a lag~1.
I have been posting slls of late but two function 'shells', CMWC( Opt OptionalLag As Long ) and CWMCSeed( Opt Text As String ), are presented here, so you can kick them around to suit your requirements.
The code presented here uses a CMWC base 2^32-1 with lag~r where r = 2, 4, 8, ..., 4096. With r = 4096 we have a period of 2^131086 compared with 2^19937 for the Mersenne Twister PRNG.
The code assumes at least Windows XP SP3 and the latest compilers but can be re-written to accomodate all versions of Windows and compilers but with more effort.
In Marsaglia's post here we have
Code:
static unsigned long Q[4096],c=362436; /* choose random initial c<809430660 and */
/* 4096 random 32-bit integers for Q[] */
unsigned long CMWC4096(void){
unsigned long long t, a=18782LL;
static unsigned long i=4095;
unsigned long x,r=0xfffffffe;
i=(i+1)&4095;
t=a*Q[i]+c;
c=(t>>32); x=t+c; if(x<c){x++;c++;}
return(Q[i]=r-x); }
The principle of CMWC is x = b - 1 - x where b is the base and x, on the right hand side, is the number generated by the MWC. With b = 2^32-1 we get x = 2^32 - 2 - x. 2^32 - 2 is 0xfffffffe; the value used for r in 'return(Q[i]=r-x'. The x's generated by MWC base 2^32-1 will be in the range [0,2^32-2]. If we use x = 0 and x = 2^32-2 in the CMWC we get [2^32-2,0]. So, the above is using a base of 2^32-1.
The term x=t+c had me fooled for a while, I was expecting a modulo.
In general, j = k*(j\k) + j MOD k ie quotient + remainder
Rearranging, j MOD k = j - k*(j\k)
With k = 2^32-1 we get
j MOD (2^32-1) = j - (2^32-1)*(j\(2^32-1))
= j - 2^32*(j\(2^32-1)) + j\(2^32-1)
Substituting Marsaglia's variables we get
x = t MOD (2^32-1) = t - 2^32*(t\(2^32-1)) + t\(2^32-1)
= t - 2^32*c + c ; with c = t\(2^32-1)
= t + c
Since 2^32*c = 0 because c is 32 bit.
I had not seen that before - it is neat isn't it?
This further confirms a base of 2^32-1.
In the following code x = t Mod &H0FFFFFFFF is used rather than x = t + c because it was found to be marginally faster; we are using Register t As Ext. I should imagine that x = t + c will come into its own with an assembler version.
So, why is c=(t>>32) being used instead of c = t\&H0FFFFFFFF? ie t\(2^32-1).
Answers on a post card to ......
In the following code c = t\&H0FFFFFFFF is used.
I should add that every implementation seen, from using C/C++, PureBASIC and Java, port Marsaglia's post as written - no questions asked.
At the moment I don't understand 'if(x<c){x++;c++;}' but have accepted it. More thought is required. Perhaps someone can put me wise on that?
For byte stream testing purposes we need to map [0,2^32-2] into [0,2^32-1] otherwise the fourth byte of a dword dump will be in the range [0,254] and not [0,255].
In addition to George Marsaglia's DIEHARD tests I use ent.exe from Fourmilab.
Here are the results from ten tests of 15MB dumps; five done one day and another five done another day.
Code:
1 2 3 4 5 6 7 8 9 10 Entropy 7.999989 7.999989 7.999989 7.999990 7.999988 7.999989 7.999988 7.999988 7.999988 7.999990 Chi square 86.67% 70.67% 83.99% 90.55% 22.72% 83.27% 55.25% 47.67% 56.03% 94.64% Arithmetic Mean 127.4922 127.4959 127.4820 127.5184 127.5042 127.5056 127.5123 127.4812 127.5121 127.5221 Monte Carlo 0.01% 0.05% 0.02% 0.04% 0.03% 0.01% 0.01% 0.02% 0.02% 0.05% SCC 0.000210 0.000120 0.000159 -0.000144 -0.000528 0.000407 -0.000089 0.000040 0.000266 0.000254
With the chi square, in particular, non-random behaviour will put us into the regions < 10% and > 90%. We only have two in the above and they would qualify as almost suspect. The more tests we do the more likely one would qualify as suspect. There is nothing unusual with that.
The conclusion here is that there is no evidence to contradict the null hypothesis that CMWC produces good quality pseudo random numbers.
Compare with these tests of 15MB dumps.
Code:
1 2 3 Entropy 7.999990 7.999989 7.999988 Chi square 97.13% 80.63% 41.68% Arithmetic Mean 127.5034 127.5028 127.46772 Monte Carlo 0.00% 0.02% 0.06% SCC 0.000156 0.000354 -0.000165
That chi square of 97.13% puts us into the suspect region but it may be a while before we see that again.
Obviously, there are many ways we can populate Q() with the initial random seeds. In the following code the CryptGenRandom API is used. It is easy to implement and, by definition, does not require seeding itself.
Best practice would call the function CMWC once in an application's initial stages so as to initialize the generator. It has an optional lag parameter which defaults to lag~256 if not given or if an invalid lag is given.
On my machine I am getting about 10 million dwords a second, irrespective of which lag is used - a bit of assembler would help.
So, which lag to use? I had reservations about a small multipler but they are not borne out in tests. There are some scenarios where exceptionally large periods are required where even Mersenne Twister is left wanting. A lag of 1024 and above exceed Mersenne Twister. Devices limited in memory would benefit from low memory Q() usage. A lag~64 needs only 256 bytes for Q() and has a period of 2^2077. You decide.
Post any comments here. I have one of my own.
Example code:
Code:
#Compile Exe
#Dim All
#Include "WIN32API.INC"
Global hProv As Dword
Function PBMain() As Long
Local i, n As Long
Local temp As Dword
Local msg As String
hProv = GetCSPHandle
If hProv = 0 Then
MessageBox 0, "Failed to get Cryptographic Service Provider handle" + $CrLf + _
"Probable cause: Operating system < XP SP3", "CMWC", %MB_IconError Or %MB_ApplModal Or %MB_Topmost
Exit Function
End If
CMWC(4096) ' Initialize generator ' takes 0.824ms On my machine
n = 20
' n = 15*1024*1024\4
' Open "15MB04096.txt" For Binary As #1
For i = 1 To n
msg = msg + Str$(CMWC) + $CrLf
' temp = CMWC()
' Put #1, , temp
Next
' Close #1
MsgBox msg
' MsgBox "Done"
If hProv Then CryptReleaseContext hProv, 0
End Function
Function CMWC( Opt OptionalLag As Long ) As Dword
Dim Q() As Static Dword ' Random x seeds
Static i, j, Multiplier, BeenHere As Long
Static cary As Dword ' Random And < Multiplier
Register x As Dword ' MWC Random number generated
Register t As Ext
If IsFalse BeenHere Then ' Initialize Q() And carry
Local Lag As Long
If IsMissing(OptionalLag) Then
Lag = 256
Else
Lag = OptionalLag
End If
i = Lag -1
j = i
ReDim Q(0 To i) As Static Dword
CryptGenRandom( hProv, 4 * Lag, VarPtr( Q(0) ) ) ' Random x seeds
Local k As Long
For k = 0 To i
If Q(k) = 4294967295 Then Decr Q(k)
' x seeds should be less than 2^32-1
Next
Select Case As Long Lag
Case 4
Multiplier = 987654366 ' 2^157
Case 8
Multiplier = 987651386 ' 2^285
Case 16
multiplier = 987765178 ' 2^541
Case 32
Multiplier = 987655670 ' 2^1053
Case 64
Multiplier = 987651206 ' 2^2077
Case 128
Multiplier = 987688302 ' 2^4125
Case 256
Multiplier = 987662290 ' 2^8221
Case 512
Multiplier = 123462658 ' 2^16410
Case 1024
Multiplier = 5555698 ' 2^32790
Case 2048
Multiplier = 1030770 ' 2^65559
Case 4096
Multiplier = 18782 ' 2^131086
Case Else
MsgBox "Invalid lag given." + $CrLf + $CrLf + "A 256 lag will be used.", %MB_Ok, "CMWC"
Multiplier = 987662290 ' 2^8221
End Select
CryptGenRandom(hProv, 4, VarPtr(cary))
cary = cary * (Multiplier-1)/4294967295 ' Make cary < Multiplier
BeenHere = %True
End If
' { Assembler here would be nice
i = ( i + 1 ) And j ' Cycle around j
t = Multiplier * Q(i) + cary
cary = t \ &H0FFFFFFFF ' Base = 2^32-1
' x = t + cary
x = t Mod &H0FFFFFFFF
' If ( x < cary ) Then ' See next post why commented out
' Incr x : Incr cary
' End If
Q(i) = &H0FFFFFFFE - x ' The CMWC aspect (Base-1) - x
' }
' Q(i) = Q(i) * 1.000000000232831## ; used when creating test files
Function = Q(i)
End Function
' ***************************************************
Function GetCSPHandle As Long
' We have a failure If the Function returns %False
Local hCSP As Dword
If CryptAcquireContext( hCSP, ByVal %Null, ByVal %Null, %PROV_RSA_AES, 0 ) = %False Then
' Unable To Find a Default For %PROV_RSA_AES so force the issue first For XP
If CryptAcquireContext( hCSP, ByVal %Null, $MS_ENH_RSA_AES_PROV_XP_A, %PROV_RSA_AES, %CRYPT_VERIFYCONTEXT ) = %False Then
' Still failed so force the issue For Vista And Win7
If CryptAcquireContext( hCSP, ByVal %Null, $MS_ENH_RSA_AES_PROV_A, %PROV_RSA_AES, %CRYPT_VERIFYCONTEXT ) <> %False Then
Function = hCSP
End If
Else
Function = hCSP
End If
Else ' found Default
Function = hCSP
End If
End Function
' ***************************************************
***********************************************************************************
There are situations where we need to have the same sequence per instance of generation along the lines of 'Randomize [number]'. The function CMWC will not do that as it stands. The CryptGenRandom aspect will need to be replaced.
One way is presented with the function CMWCSeed along the lines of 'Randomize [alphanumeric]' which, of course, we do not have with PB.
The principle is to take an alphanumeric seed and hash it with, say, SHA256 to populate a lag~8 Q(). This has a period of 2^285. To put that into perspective 2^285/2^32 = 2^253 ~ 10^76 times the period of RND. A SHA512 approach will give ~ 10^162 times the period of RND.
The alphanumeric parameter is optional. If not given "The time has come the walrus said...." is used; and why not?
Obviously, there are many ways to populate Q() with a seeded approach - using a hash is an easy way.
Example code:
Code:
#Compile Exe
#Dim All
#Include "WIN32API.INC"
Global hProv As Dword
Function PBMain() As Long
Local i, n As Long
Local msg As String
hProv = GetCSPHandle
If hProv = 0 Then
MessageBox 0, "Failed to get Cryptographic Service Provider handle" + $CrLf + _
"Probable cause: Operating system < XP SP3", "CMWC", %MB_IconError Or %MB_ApplModal Or %MB_Topmost
Exit Function
End If
CWMCSeed("PowerBASIC Compile without compromise") ' Takes 0.043ms On my machine
n = 20
For i = 1 To n
msg = msg + Str$(CWMCSeed) + $CrLf
Next
MsgBox msg
If hProv Then CryptReleaseContext hProv, 0
End Function
Function CWMCSeed( Opt Text As String ) As Dword
Dim Q() As Static Dword ' Random x seeds
Static i, BeenHere As Long
Static cary As Dword
Register x As Dword ' MWC Random number generated
Register t As Ext
If IsFalse BeenHere Then ' Initialize Q()
Local TextToCrunch As String
If IsMissing(Text) Then
TextToCrunch = "The time has come the walrus said...."
Else
TextToCrunch = Text
End If
Local hHash, hashSize As Long
i = 7
ReDim Q(0 To i) As Static Dword
CryptCreateHash( hProv, %CALG_SHA_256, 0, 0, hHash )
CryptHashData( hHash, ByVal StrPtr(TextToCrunch), Len(TextToCrunch), 0 )
CryptGetHashParam( hHash, %HP_HASHSIZE, VarPtr(hashSize), 4, 0 )
CryptGetHashParam( hHash, %HP_HASHVAL, VarPtr( Q(0) ), hashSize, 0 )
Local k As Long
For k = 0 To i
If Q(k) = 4294967295 Then Decr Q(k)
' x seeds should be less than 2^32-1
Next
cary = 123456789 ' Carved In stone but could be anything As
' Long As it Is Fixed And less than 9877651386
BeenHere = %True
End If
i = ( i + 1 ) And 7
t = 987651386 * Q(i) + cary
cary = t \ &H0FFFFFFFF ' Base = 2^32-1
'x = t + cary
x = t Mod &H0FFFFFFFF
' If ( x < cary ) Then ' See next post why commented out
' Incr x : Incr cary
' End If
Q(i) = &H0FFFFFFFE - x ' The CMWC aspect (Base-1) - x
Function = Q(i)
End Function
' *****************************************************************************************************
Function GetCSPHandle As Long
' We have a failure If the Function returns %False
Local hCSP As Dword
If CryptAcquireContext( hCSP, ByVal %Null, ByVal %Null, %PROV_RSA_AES, 0 ) = %False Then
' Unable To Find a Default For %PROV_RSA_AES so force the issue first For XP
If CryptAcquireContext( hCSP, ByVal %Null, $MS_ENH_RSA_AES_PROV_XP_A, %PROV_RSA_AES, %CRYPT_VERIFYCONTEXT ) = %False Then
' Still failed so force the issue For Vista And Win7
If CryptAcquireContext( hCSP, ByVal %Null, $MS_ENH_RSA_AES_PROV_A, %PROV_RSA_AES, %CRYPT_VERIFYCONTEXT ) <> %False Then
Function = hCSP
End If
Else
Function = hCSP
End If
Else ' found Default
Function = hCSP
End If
End Function
' *****************************************************************************************************
Comment