FreeBASIC's PRNG #2

General FreeBASIC programming questions.
Post Reply
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: FreeBASIC's PRNG #2

Post by deltarho[1859] »

@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.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: FreeBASIC's PRNG #2

Post by deltarho[1859] »

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!
paul doe
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: FreeBASIC's PRNG #2

Post by paul doe »

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.
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:Well, well, well - what do you know? 1TB and only three lowest ranking anomalies.
Ha, looks can be deceptive, no? =D
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.
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: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'.
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: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.
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 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.
Sure, why not? Meet the Middle Square Weyl Sequence RNG, Tinkerer's Version:

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()
Let me know if there's a snafu somewhere, I shall fix it immediately. Now, have fun with it! =D
deltarho[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!
Glad to know that you're getting your rush with this cute little monster =D
Last edited by paul doe on Sep 12, 2018 23:52, edited 1 time in total.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: FreeBASIC's PRNG #2

Post by deltarho[1859] »

paul doe wrote:Let me know if there's a snafu somewhere, I shall fix it immediately. Now, have fun with it! =D
I will have a look tomorrow. A quick glance and your latest code's initialization is similar to PCG32II.

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.
Glad to know that you're getting your rush with this cute little monster =D
Many will be rolling their eyes at my reaction but each to their own. <smile>
paul doe
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: FreeBASIC's PRNG #2

Post by paul doe »

deltarho[1859] wrote:I will have a look tomorrow. A quick glance and your latest code's initialization is similar to PCG32II.
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: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.
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 =D
deltarho[1859] wrote:Many will be rolling their eyes at my reaction but each to their own. <smile>
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 =D
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: FreeBASIC's PRNG #2

Post by dodicat »

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.

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

  
My chi squared distance is 203.78
MrSwiss
Posts: 3910
Joined: Jun 02, 2013 9:27
Location: Switzerland

Re: FreeBASIC's PRNG #2

Post by MrSwiss »

paul doe wrote:Let me know if there's a snafu somewhere, I shall fix it immediately. Now, have fun with it! =D
I think there is one:
in the paper it states that, last bit must be set always ...
You are doing: And lsb ... which should be, IMO:
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.
paul doe
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: FreeBASIC's PRNG #2

Post by paul doe »

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.
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?
Last edited by paul doe on Sep 12, 2018 22:50, edited 1 time in total.
paul doe
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: FreeBASIC's PRNG #2

Post by paul doe »

MrSwiss wrote:I think there is one:
in the paper it states that, last bit must be set always ...
You are doing: And lsb ... should be, IMO:
Or lsb (to guarantee the requirement, as outlined below), NO?
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:

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
Naturally, you can also do it as you suggest (which is also the method suggested by the author)
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: FreeBASIC's PRNG #2

Post by deltarho[1859] »

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?
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.

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.
paul doe
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: FreeBASIC's PRNG #2

Post by paul doe »

deltarho[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.
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 =D
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: FreeBASIC's PRNG #2

Post by deltarho[1859] »

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

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
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: FreeBASIC's PRNG #2

Post by deltarho[1859] »

@paul doe

Testing your code.

I got an immediate PractRand FAILure when I used:

Code: Select all

var randomNumber = MsWs(987654, 0, 8)
Of course, the given sequence is not vaild and isValidSequence(8) told me so.

So, I tried

Code: Select all

var randomNumber = MsWs(987654, 0, 8ull shl 32 or 1)
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.

Code: Select all

var randomNumber = MsWs(987654, 1, 8ull shl 32 or 1) 
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>
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: FreeBASIC's PRNG #2

Post by jj2007 »

deltarho[1859] wrote:From CryptoRNDII
I've tried to examine what this code does; IMHO it makes the result less random than before:

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
Output:

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
The punpckldq is ok, it just creates a 64-bit random number from 2x32 bits (but movlps xmm0, qword ptr [eax] is 10 bytes shorter and does exactly the same). But the rest? The 3FF is a constant, for example...
dafhi
Posts: 1640
Joined: Jun 04, 2005 9:51

Re: FreeBASIC's PRNG #2

Post by dafhi »

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 )
Last edited by dafhi on Sep 13, 2018 19:10, edited 1 time in total.
Post Reply