FreeBASIC's PRNG #2
-
- Posts: 4292
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
@paul doe
Paul, your understanding is slightly off.
m_s is not the seed but the sequence or stream number.
m_x is the seed.
m_w is the Weyl sequence.
The last random number is stored in m_x. However, m_w is a stepping sequence and the value of m_s is added to it on each iteration.
So, we can regard m_x and m_w as the state vector.
We have a 32 bit generator with a period of 2^64 and a 4*32-bit state vector. By comparison Mersenne Twister has a large state vector of 624*32-bits; considered by some as being a tad 'heavy'.
Your code then is creating a random sequence number and uses zero for both the seed and initial Weyl sequence. The initial x * x will then be zero but we avoid the 'zero mechanism', mentioned in the paper, by adding the Weyl sequence. Your initial Weyl sequence is zero but this is not a problem because it is 'stepped up' by the sequence number before we start to do any shifting.
Your code then should allow either a fixed sequence value or a random sequence value, via m_s, and a seed, m_x, which may be zero. The initial Weyl sequence, m_w, is an unsigned 64-bit integer and may be zero. The Weyl sequence is analogous to the Initialization Vector (IV) used in AES.
Paul, your understanding is slightly off.
m_s is not the seed but the sequence or stream number.
m_x is the seed.
m_w is the Weyl sequence.
The last random number is stored in m_x. However, m_w is a stepping sequence and the value of m_s is added to it on each iteration.
So, we can regard m_x and m_w as the state vector.
We have a 32 bit generator with a period of 2^64 and a 4*32-bit state vector. By comparison Mersenne Twister has a large state vector of 624*32-bits; considered by some as being a tad 'heavy'.
Your code then is creating a random sequence number and uses zero for both the seed and initial Weyl sequence. The initial x * x will then be zero but we avoid the 'zero mechanism', mentioned in the paper, by adding the Weyl sequence. Your initial Weyl sequence is zero but this is not a problem because it is 'stepped up' by the sequence number before we start to do any shifting.
Your code then should allow either a fixed sequence value or a random sequence value, via m_s, and a seed, m_x, which may be zero. The initial Weyl sequence, m_w, is an unsigned 64-bit integer and may be zero. The Weyl sequence is analogous to the Initialization Vector (IV) used in AES.
-
- Posts: 4292
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
BREAKING NEWS
I left PractRand running and we got to 2TB without any further issues. I pulled out to give my CPU a well-earned rest. <smile>
WOW!
I left PractRand running and we got to 2TB without any further issues. I pulled out to give my CPU a well-earned rest. <smile>
WOW!
Re: FreeBASIC's PRNG #2
Yeah, crazy stuff huh? I found it interesting, since it looks quite straightforward (and even naïve to me, but what do I know =D)deltarho[1859] wrote:I read the paper, looked at your 'function one() as ulong' and figured Widynski was having a laugh especially since von Neumann died in 1957.
Ha, looks can be deceptive, no? =Ddeltarho[1859] wrote:Well, well, well - what do you know? 1TB and only three lowest ranking anomalies.
Ok, understood. The confusion came because of the usage of the word 'seed' in the paper by the author, to refer to the 'sequence', and the fact that there's a 'seeds.h' file in the package I linked before which contains 25000 valid such sequences. Much clearer now, thanks.deltarho[1859] wrote:Paul, your understanding is slightly off.
m_s is not the seed but the sequence or stream number.
m_x is the seed.
m_w is the Weyl sequence.
The last random number is stored in m_x. However, m_w is a stepping sequence and the value of m_s is added to it on each iteration.
So, we can regard m_x and m_w as the state vector.
Yeah, the thing is tiny. I liked it precisely for this characteristic, because you can then instantiate as many of them as you want, and the memory footprint is minimal. This would allow me to use it to great effect in, say, procedural content generation (my main usage for such generators). Having a procedurally generated object to have its own RNG embedded in it is a powerful concept.deltarho[1859] wrote:We have a 32 bit generator with a period of 2^64 and a 4*32-bit state vector. By comparison Mersenne Twister has a large state vector of 624*32-bits; considered by some as being a tad 'heavy'.
Yes, I just ported the code 'as is'. I think that the author's intent was precisely that: to show that very minimal state is needed to generate good numbers, and it also doesn't need any kind of warm up, as far as I can tell.deltarho[1859] wrote:Your code then is creating a random sequence number and uses zero for both the seed and initial Weyl sequence. The initial x * x will then be zero but we avoid the 'zero mechanism', mentioned in the paper, by adding the Weyl sequence. Your initial Weyl sequence is zero but this is not a problem because it is 'stepped up' by the sequence number before we start to do any shifting.
Sure, why not? Meet the Middle Square Weyl Sequence RNG, Tinkerer's Version:deltarho[1859] wrote:Your code then should allow either a fixed sequence value or a random sequence value, via m_s, and a seed, m_x, which may be zero. The initial Weyl sequence, m_w, is an unsigned 64-bit integer and may be zero. The Weyl sequence is analogous to the Initialization Vector (IV) used in AES.
Code: Select all
/'
The original code can be found at:
https://mswsrng.wixsite.com/rand
uint64_t x = 0, w = 0, s = 0xb5ad4eceda1ce2a9;
inline static uint32_t msws() {
x *= x;
x += (w += s);
return x = (x>>32) | (x<<32);
}
'/
#define hiqword( x ) ( ( cast( ulongint, x ) and &hffffffff00000000 ) shr 32 )
#define loqword( x ) ( cast( ulongint, x ) and &h00000000ffffffff )
#define lsb &b1 '' Least-significant bit
type MsWs
public:
declare constructor()
declare constructor( _
byval as ulongint )
declare constructor( _
byval as ulongint, _
byval as ulongint )
declare constructor( _
byval as ulongint, _
byval as ulongint, _
byval as ulongint )
declare destructor()
declare property initialSeed() as ulongint
declare property initialSequence() as ulongint
declare property initialWeilSequence() as ulongint
declare property seed() as ulongint
declare property sequence() as ulongint
declare property weilSequence() as ulongint
declare function one() as ulong
declare static function validRandomSequence() as ulongint
declare static function isValidSequence( _
byval as ulongint ) as boolean
private:
'' Just for testing purposes
m_initialSeed as ulongint
m_initialWeilSequence as ulongint
m_initialSequence as ulongint
'' The real state of the PRNG
m_x as ulongint
m_weilSeq as ulongint
m_seq as ulongint
end type
constructor MsWs()
this.constructor( 0, 0, 0 )
end constructor
constructor MsWs( _
byval aSeed as ulongint )
this.constructor( aSeed, 0, 0 )
end constructor
constructor MsWs( _
byval aSeed as ulongint, _
byval aWeilSeq as ulongint )
this.constructor( aSeed, aWeilSeq, 0 )
end constructor
constructor MsWs( _
byval aSeed as ulongint, _
byval aWeilSeq as ulongint, _
byval aSeq as ulongint )
m_x = aSeed
m_weilSeq = aWeilSeq
if( aSeq = 0 ) then
m_seq = validRandomSequence()
else
m_seq = aSeq
end if
m_initialSeed = m_x
m_initialWeilSequence = m_weilSeq
m_initialSequence = m_seq
end constructor
destructor MsWs()
end destructor
property MsWs.initialSeed() as ulongint
return( m_initialSeed )
end property
property MsWs.initialWeilSequence() as ulongint
return( m_initialWeilSequence )
end property
property MsWs.initialSequence() as ulongint
return( m_initialSequence )
end property
property MsWs.seed() as ulongint
return( m_x )
end property
property MsWs.weilSequence as ulongint
return( m_weilSeq )
end property
property MsWs.sequence() as ulongint
return( m_seq )
end property
function MsWs.isValidSequence( _
byval aSequence as ulongint ) as boolean
return( cbool( _
( hiqword( aSequence ) <> 0 ) andAlso _
( loqword( aSequence ) and lsb ) ) )
end function
function MsWs.validRandomSequence() as ulongint
/'
According to Bernard Widynski (the creator of this PRNG), the sequence
must be a number that has the upper 32-bits non-zero and the least
significant bit set to 1. There are 2 ^ 63 numbers that satisfy this
property, and 25000 such numbers are given in the C implementation.
Here, I just return a random sequence that satisfies that property.
'/
randomize( , 5 )
dim as ulongint aSequence
do while( not isValidSequence( aSequence ) )
aSequence = ( culngint( rnd() * culngint( -1 ) ) shr 32 ) or _
culngint( rnd() * culngint( -1 ) )
loop
return( aSequence )
end function
function MsWs.one() as ulong
m_x *= m_x : m_weilSeq += m_seq : m_x += m_weilSeq
m_x = ( m_x shr 32 ) or ( m_x shl 32 )
return( m_x )
end function
/'
Unit test
var randomNumber = MsWs( 0, 0, &hb5ad4eceda1ce2a9 ) '' <-- To test equivalence with C code
var randomNumber = MsWs() '' <-- The default setting creates a random valid sequence
var randomNumber = MsWs( aSeed ) '' <-- Using a seed
var randomNumber = MsWs( aSeed, aWeilSequence ) '' <-- Using a seed and a custom Weil Sequence
var randomNumber = MsWs( aSeed, aWeilSequence, aSequence ) '' <-- The three parameters specified
You can omit a parameter simply by passing 0
var randomNumber = MsWs( aSeed, 0, aSequence )
var randomNumber = MsWs( 0, aWeilSequence, aSequence )
Or ask the MsWs class to generate a valid sequence for you:
var randomNumber = MsWs( 0, 0, MsWs.validRandomSequence() )
You can also ask it if a sequence is valid (and act accordingly otherwise):
dim as ulongint aSequence = iif( _
MsWs.isValidSequence( &Hdeadf00d ), _ '' <-- is this one valid?
&Hdeadf00d, MsWs.validRandomSequence() )
To clone an instance, create another one, like this:
var anotherRandomNumber = MsWs( _
randomNumber.seed, _
randomNumber.weilSequence, _
randomNumber.sequence )
There. Now you have 2 exact copies of the same RNG, and should return identical
results when requesting one() random number from them. Useful for testing them in
parallel (so you can test, say, 4 or 5 parameter sets at once)
Of course, you can also clone it with any of the three parameters:
var anotherRandomNumber = MsWs( _ '' <-- Same seed and Weil Sequence, different sequence
randomNumber.seed, _
randomNumber.weilSequence )
var anotherRandomNumber = MsWs( _ '' <-- You get the idea
someOtherSeed, _
randomNumber.weilSequence _
randomNumber.sequence )
'/
var randomNumber = MsWs()
for i as integer = 1 to 100
? randomNumber.one()
next
?
? "How to clone your instance..."
?
var anotherRandomNumber = MsWs( _
randomNumber.seed, _
randomNumber.weilSequence, _
randomNumber.sequence )
for i as integer = 1 to 100
? randomNumber.one(), anotherRandomNumber.one()
next
sleep()
Glad to know that you're getting your rush with this cute little monster =Ddeltarho[1859] wrote:BREAKING NEWS
I left PractRand running and we got to 2TB without any further issues. I pulled out to give my CPU a well-earned rest. <smile>
WOW!
Last edited by paul doe on Sep 12, 2018 23:52, edited 1 time in total.
-
- Posts: 4292
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
I will have a look tomorrow. A quick glance and your latest code's initialization is similar to PCG32II.paul doe wrote:Let me know if there's a snafu somewhere, I shall fix it immediately. Now, have fun with it! =D
Once the 'sea trials' have been completed I shall add single(extended), double and range to start with so that we can 'drop into' our FB code in place of RND.
Many will be rolling their eyes at my reaction but each to their own. <smile>Glad to know that you're getting your rush with this cute little monster =D
Re: FreeBASIC's PRNG #2
Indeed. O'Neil mentions Weyl sequences in her blog (which is quite a good read, BTW) and how they're used in a similar fashion in the PCG family of generators. Concretely, she states that 'The extension array is a funky counter akin to a Weyl sequence (each array element is like a digit of a counter)', so that both algorithms share commonalities should not come as a surprise.deltarho[1859] wrote:I will have a look tomorrow. A quick glance and your latest code's initialization is similar to PCG32II.
Yes, of course. This is the 'test setup', so to speak. Once all is good, we can lay it out differently (a class for me, some functions for other folks who prefer the procedural approach). The code is trivial in any case =Ddeltarho[1859] wrote:Once the 'sea trials' have been completed I shall add single(extended), double and range to start with so that we can 'drop into' our FB code in place of RND.
Yep, but a good, fast and reliable RNG is crucial for anybody. Making pretty pictures appear on screen seems like a vanity affair compared with this stuff =Ddeltarho[1859] wrote:Many will be rolling their eyes at my reaction but each to their own. <smile>
Re: FreeBASIC's PRNG #2
For curiosity I ran my own random string generator through deltarho's chisquare2g loop.
It is via a string, so it takes 4 or 5 seconds.
My chi squared distance is 203.78
It is via a string, so it takes 4 or 5 seconds.
Code: Select all
#Include "file.bi"
Function rndX overload (s1 As String) As String
#macro GetNumber
#define range(f,l) Int(Rnd*((l+1)-(f))+(f))
s[0]=range(48,s1[0])
For n As Long = 1 To L-1
s[n]=range(48,57)
Next
#endmacro
#macro compare(n1,n2,ans)
lenn1=Len(n1):lenn2=Len(n2)
If lenn1 > lenn2 Then ans=-1:Goto lbl
If lenn1 < lenn2 Then ans=0:Goto lbl
If n1 > n2 Then ans = -1 Else ans= 0
lbl:
#endmacro
Dim As Long L=Len(s1),ans=1
dim as long lenn1,lenn2
Dim As String s=String(L,0)
While ans
GetNumber
compare(s,s1,ans)
Wend
Return Ltrim(s,"0")
End Function
If FileExists( "16MBs.txt" ) Then Kill "16MBs.txt"
Dim As Ulong i, n = 4*1024*1024
dim as double t=timer
Open "16MBs.txt" For Binary As #1
For i = 1 To n
Put #1,, Cast( Ulong, valulng(rndX(str(2^32)) ))
Next
Close #1
print "done",timer-t
sleep
Re: FreeBASIC's PRNG #2
I think there is one:paul doe wrote:Let me know if there's a snafu somewhere, I shall fix it immediately. Now, have fun with it! =D
in the paper it states that, last bit must be set always ...
You are doing:
Or lsb (to guarantee the requirement, as outlined below), NO?
When generating seeds for parallel usage, one might simply randomize
the s value with a 64-bit random number, verify that the upper 32 bits
are non-zero, and set the least significant bit to 1. This will with great
probability produce all different seeds. If it is necessary to guarantee that
all are different, one might pre-compute a set of random s values and verify
that all were different and then use this set when creating streams.
Last edited by MrSwiss on Sep 12, 2018 22:37, edited 1 time in total.
Re: FreeBASIC's PRNG #2
I would like to modify this fella to return 64-bit numbers also. Seems straightforward enough, but I'll need your feedback, since I'm not up to speed on the topic. What do you think? Is it feasible/needed?deltarho[1859] wrote:Once the 'sea trials' have been completed I shall add single(extended), double and range to start with so that we can 'drop into' our FB code in place of RND.
Last edited by paul doe on Sep 12, 2018 22:50, edited 1 time in total.
Re: FreeBASIC's PRNG #2
Indeed, the paper specifies that. However, if you look closely, it checks (not sets) the lsb and, if not set, simply generates another random number until both properties are satisfied. Thus, the requirement is guaranteed. You can check the correctness of the code simply by passing invalid input to the 'isValidSequence()' method, since it's static:MrSwiss wrote:I think there is one:
in the paper it states that, last bit must be set always ...
You are doing:Andlsb ... should be, IMO:
Or lsb (to guarantee the requirement, as outlined below), NO?
Code: Select all
? MsWs.isValidSequence( &h00000000000b1eed ) '' Returns 'false', upper 32-bits are zero
? MsWs.isValidSequence( &hf00d000000000000 ) '' Returns 'false', lsb is zero
? MsWs.isValidSequence( &hb5ad4eceda1ce2a9 ) '' Returns 'true'; this is the sequence provided in the C code
-
- Posts: 4292
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
Feasible?: I would say so and 'off the top of my head', I reckon that we would need either direct 64-bit arithmetic or emulation.paul doe wrote:I would like to modify this fella to return 64-bit numbers also. Seems straightforward enough, but I'll need your feedback, since I'm not up to speed on the topic. What do you think? Is it feasible/needed?
Needed?: I am not sure on that one. A lot of academics get involved in solutions in search of a problem.
Manufacturer: I have a problem.
Professor: Explain it to me.
Manufacturer: Blah, blah, blah.
Professor: Turn around, the solutions shelf is the third shelf up and you need the one fifth from the right. Do you want a plastic bag for it?
Manufacturer: No, I have brought my own.
Professor: Oh, well done.
Re: FreeBASIC's PRNG #2
I was thinking of this: proceed as normal but with two independent seeds and sequences. Then, after each one has been mixed, we combine them into a 64-bit number. Do you think this crude approach is workable? We don't have native 128-bit operations in FB yet, so we have to make do with what we have. However, if I'm sounding too noobish, please, slap me in the forehead, because indeed I am =Ddeltarho[1859] wrote:Feasible?: I would say so and 'off the top of my head', I reckon that we would need either direct 64-bit arithmetic or emulation.
-
- Posts: 4292
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
I have recently added randD to PCG32II. This takes two 32-bit random numbers and converts them into a Double resulting in the Double having 53-bit granularity. I have been doing this with CryptoRNDII from its inception.
Now some may argue that since PRNGs of the form xn+1 = f(xn) then there will exist a correlation between xn+1 and xn. That is true but a top-drawer generator will see that correlation as almost non-existent.
In the ENT test that I did with the MSWS RNG the fourth metric, serial correlation, was -0.000097. Most top drawer generators see at least three zeros after the decimal point.
Admittedly if we took two 32-bit random numbers from two independent generators then the serial correlation coefficient should be zero. I take the view that with a figure as small as -0.000097 then two 32-bit random numbers from the same generator will be fine.
Either way, my response to your proposal is: "Yes, it is workable"
There is no need to 'mix' them. If two numbers are random then their concatenation will be random. Mixing them will 'cost money'.
From CryptoRNDII
Now some may argue that since PRNGs of the form xn+1 = f(xn) then there will exist a correlation between xn+1 and xn. That is true but a top-drawer generator will see that correlation as almost non-existent.
In the ENT test that I did with the MSWS RNG the fourth metric, serial correlation, was -0.000097. Most top drawer generators see at least three zeros after the decimal point.
Admittedly if we took two 32-bit random numbers from two independent generators then the serial correlation coefficient should be zero. I take the view that with a figure as small as -0.000097 then two 32-bit random numbers from the same generator will be fine.
Either way, my response to your proposal is: "Yes, it is workable"
There is no need to 'mix' them. If two numbers are random then their concatenation will be random. Mixing them will 'cost money'.
From CryptoRNDII
Code: Select all
Private Function CryptoD As Double ' [0,1)
If ptrBuffer >= SwitchBufferCriteria Then
SwitchBuffer
End If
Asm
mov eax, dword Ptr [ptrBuffer]
movd xmm0, [eax] ' first 32-bit random number
movd xmm1, [eax + 4] ' second 32-bit random number
punpckldq xmm0, xmm1
psrlq xmm0, 12
mov eax, 1
cvtsi2sd xmm1, eax
por xmm0, xmm1
subsd xmm0, xmm1
movq [Function], xmm0
End Asm
ptrBuffer += 8
End Function
-
- Posts: 4292
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
@paul doe
Testing your code.
I got an immediate PractRand FAILure when I used:
Of course, the given sequence is not vaild and isValidSequence(8) told me so.
So, I tried
after isValidSequence(8ull shl 32 or 1) told me the sequence was OK.
However, PractRand still FAILed immediately so it looks like passing zero for the weilSequnce is causing a problem.
also FAILed so the problem is with the second parameter.
I am going through the code but you will probably spot the problem before me, being the author of the code.
A newbie would greatly benefit from reading your code as it is a master class, no pun intended, on using Constructor and Property. Yours truly is benefitting. <smile>
Testing your code.
I got an immediate PractRand FAILure when I used:
Code: Select all
var randomNumber = MsWs(987654, 0, 8)
So, I tried
Code: Select all
var randomNumber = MsWs(987654, 0, 8ull shl 32 or 1)
However, PractRand still FAILed immediately so it looks like passing zero for the weilSequnce is causing a problem.
Code: Select all
var randomNumber = MsWs(987654, 1, 8ull shl 32 or 1)
I am going through the code but you will probably spot the problem before me, being the author of the code.
A newbie would greatly benefit from reading your code as it is a master class, no pun intended, on using Constructor and Property. Yours truly is benefitting. <smile>
Re: FreeBASIC's PRNG #2
I've tried to examine what this code does; IMHO it makes the result less random than before:deltarho[1859] wrote:From CryptoRNDII
Code: Select all
mov eax, offset ptrBuffer
movd xmm0, dword ptr [eax] ; ' first 32-bit random number
movd xmm1, dword ptr [eax + 4] ; ' second 32-bit random number
deb 4, "before", x:xmm0, x:xmm1 ; x=display as hex
punpckldq xmm0, xmm1
deb 4, "after punpckldq", x:xmm0
psrlq xmm0, 12
deb 4, "after psrlq", x:xmm0
mov eax, 1
cvtsi2sd xmm1, eax
deb 4, "after cvtsi2sd", x:xmm1
por xmm0, xmm1
deb 4, "after por", x:xmm0
subsd xmm0, xmm1
deb 4, "after subsd", x:xmm0
Code: Select all
before
x:xmm0 00000000 00000000 00000000 12345678
x:xmm1 00000000 00000000 00000000 90ABCDEF
after punpckldq x:xmm0 00000000 00000000 90ABCDEF 12345678
after psrlq x:xmm0 00000000 00000000 00090ABC DEF12345
after cvtsi2sd x:xmm1 00000000 00000000 3FF00000 00000000
after por x:xmm0 00000000 00000000 3FF90ABC DEF12345
after subsd x:xmm0 00000000 00000000 3FE21579 BDE2468A
Re: FreeBASIC's PRNG #2
I'm not sure how my randomness is, but my 'csg ii' period surpasses MsWs and PCG.
viewtopic.php?f=8&t=27007#p251801
there are ways to extend the period csg-ii, like slow the progression of 'a += 1' with
+ ( (state shr (state.cbits-1)) = 0 )
viewtopic.php?f=8&t=27007#p251801
there are ways to extend the period csg-ii, like slow the progression of 'a += 1' with
+ ( (state shr (state.cbits-1)) = 0 )
Last edited by dafhi on Sep 13, 2018 19:10, edited 1 time in total.