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