Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

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

Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

Post by deltarho[1859] »

Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs), found by Paul Doe, is doing so well in speed tests and bearing in mind that it has passed PractRand to 1TB I thought I would develop it beyond a straight Function 'drop in' to one which is thread-safe to complement PCG32II; the only thread-safe generator that I have.

The MsWs drop in Function has an issue. This was created when we 'hard-wired' s; where s should be non-zero in the upper 32 bits and 1 in the least significant bit to add to the stepping sequence w += s.
One might think that non-overlapping data would always be obtained if the generator was started with different x values but the same w and s values. It is possible for two different x values to have the same result mod 2^64 after squaring. Should this happen, exactly the same data would be produced even though the initial x values were different. To guarantee distinct data through the period of 2^64, one should use a different s value for each run. This will produce different data and prevent any occurrence of overlapping data.
The author of MsWs provides 25000 values of s; 2^14.60964047 values.

Note that we are referring to the same w and s. We are using a 'hard-wired' s so all we need do is ensure that we do not use the same w and avoid trouble that way. However, if we use a random w it is possible to generate the same w. A random w is got by using the function Get64Bit. The chance of a 'collision' is then 2^64 to 1 against which is considerably larger than 2^14.60964047. The odds of winning the US Powerball lottery is 292,201,338 to 1 against; approximately 2^28.12238754. 2^64 is greater than 2^35 x 2^28.12238754 so 2^64 is treated as 'it ain't going to happen'. <smile>

With the function MyRandomize the option MyRandomize( 0, <fixed w> ) would generate a random seed and that plus <fixed w> will lead to overlapping data. If you use MyRandomize( 0, <fixed w> ) you will not get it but will get, instead, MyRandomize(). We can use <fixed w> but only if we use <fixed seed> as well when we require to repeat a sequence.

MyRandomiize() is invoked via a Constructor so if you require full random 'seeding' then there is no need to employ MyRandomize. MyRandomize includes a small warmup. On my machine, MyRandomize takes 3.4ms to complete.

The following code has nothing like the features of PCG32II, which has been worked upon for some time, and currently only two generator functions are available namely Rand ( Ulongint ) and RandSE.

However, MsWs has more to it than meets the eye - the giveaway is in the last paragraph. With a Ulongint output, only one random number is required to give 53-bit granularity for [0,1) coming in, using badidea's loop overhead filter, at 433MHz. PCG32II, on the other, requires two Ulongs coming in at 221MHz. CryptoRNDII requires two Ulongs coming in at 442MHz. Knuth64, MsWs nearest competitor, comes in at 487MHz. Knuth64 fails PractRand at 8KB which does not compare favorably with MsWs 1TB success. In 'real life' situations the difference between 433MHz and 487MHz will not be remarkable. All of those figures are based upon 9 runs each and the median taken for FBC 1.05/gcc 5.2 64-bit. If thread safety is an issue then MsWs is 'top dog'.

I have just run the test code and got an average of 0.4999956925552035. An accuracy of five significant figures is not unusual for MsWs.

MsWs.bas ( Updated to include Range functions: 29 Oct 2018 3rd revision; 30 Oct 2018 4th revision, see MsWsParams.Range; 01 Nov 2018 5th revision, corrected a blunder in MsWsParams.Range)

Code: Select all

#define Get64Bit Cast( Ulongint, Rnd*2^64 )

Type MsWsParams
  Declare Constructor
  Declare Sub MyRandomize( Byval seed As Ulongint = 0, Byval sequence As Ulongint = 0 )
  Declare Function Rand() As Ulongint
  Declare Function RandSE() As Double
  Declare Function Range Overload ( First as Double, Last as Double ) as Double
  Declare Function Range Overload ( First as Longint, Last as Longint ) as Longint
  Seed As Ulongint
  Sequence As Ulongint
End Type

Sub MsWsParams.MyRandomize( Byval seed As Ulongint = 0, Byval sequence As Ulongint = 0 )
  Randomize , 5
  If seed = 0 Then
    This.Seed = Get64Bit
    This.Sequence = Get64Bit
  Else
    If sequence = 0 Then
      This.Seed = seed
      This.Sequence = Get64Bit
    Else ' For a sequence repeat
      This.Seed = seed
      This.Sequence = sequence
    End If
  End If
  ' Warm up generator - essential for any PRNG
  For i As Ulong = 1 To 1000
    This.rand()
  Next
End Sub

' 18446744073709551557 is the largest prime number < 2^64

Function MsWsParams.Rand() As Ulongint
  This.Seed *= This.Seed : This.Sequence += 18446744073709551557 : This.Seed += This.Sequence
  This.Seed = ( This.Seed Shr 32 ) Or ( This.Seed Shl 32 )
  Return( This.Seed )
End Function

Function MsWsParams.RandSE() As Double
  This.Seed *= This.Seed : This.Sequence += 18446744073709551557 : This.Seed += This.Sequence
  This.Seed = ( This.Seed Shr 32 ) Or ( This.Seed Shl 32 )
  Return This.Seed/2^64
End Function

Function MsWsParams.Range( First as Double, Last as Double ) as Double
  This.Seed *= This.Seed : This.Sequence += 18446744073709551557 : This.Seed += This.Sequence
  This.Seed = ( This.Seed Shr 32 ) Or ( This.Seed Shl 32 )
  Function = This.Seed/2^64*( Last - First ) + First 
End Function

' Following function from code by dafhi
Function MsWsParams.Range(First As Longint, Last As Longint) As LongInt
  Function = This.Seed/2^64 * (Last - First + 1) - .5 + First
End Function

Constructor MsWsParams
  This.MyRandomize()
End Constructor
Last edited by deltarho[1859] on Nov 03, 2018 9:30, edited 12 times in total.
dafhi
Posts: 1641
Joined: Jun 04, 2005 9:51

Re: Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

Post by dafhi »

thanks for sharing the framework

Knuth LCG experiment

Code: Select all

  state = mul * (state shr 9) + add
my "CSG_ii"

Code: Select all

    a += 1 - (state = 0)
  state = (mul * (state shl 52 or state shr 12) + add) xor a
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

Post by deltarho[1859] »

@dafhi

You should not do this:

Code: Select all

state = mul * (state shr 9) + add
You are destroying the properties of the generator. mul and add are not any old numbers and you are committing what is known as LCV (Linear Congruential Vandalism). <smile>

There is nothing wrong with conditioning to improve the output so you should do something like this:

Code: Select all

oldstate = state
state = mul * state + add ' honours generator
output = mul * (oldstate shr 9) + add ' conditioning
May I extend my heartiest congratulations. Knuth64 fails PractRand at 8KB. With my edit, your experiment got to 256KB and failed.

However, your LCV managed to get to 1MB and failed. That does not mean you were correct. You were very lucky - it could have gone the other way. You would need to go 1,073,741,824 times further to get a PractRand gold medal. I have to wear half of mine on one breast and the other half on the other breast otherwise I'd fall over. <cough>
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

Post by jj2007 »

deltarho[1859] wrote:... a PractRand gold medal. I have to wear half of mine on one breast and the other half on the other breast otherwise I'd fall over. <cough>
You made my day, thanks ;-)
dafhi
Posts: 1641
Joined: Jun 04, 2005 9:51

Re: Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

Post by dafhi »

However, your LCV managed to get to 1MB and failed
have a look at a plain jane LCG's low bits

Code: Select all

function CSG as ulong
    
  const mul = 6364136223846793005ull '' Knuth's
  const add = 1442695040888963407ull
  
  static as ulongint state = 1
  
  state = mul * state shr 0 + add
 
  return (state and 4294967295)

End function


const             COLOR_OFF = rgb(0,0,255)
const             COLOR_ON  = rgb(230,230,220)

sub PixelBits(valu as ulongint, x as short, y as short, size as short=4, bits as byte = 32)
    var lenm = size - 1
    for shift as long = bits-1 to 0 step -1
      dim as ulong col = IIF( (valu shr shift)and 1, COLOR_ON, COLOR_OFF )
      line (x,y)-(x+lenm, y+lenm),col, bf
      x += size
    Next
End Sub


sub Main

    var           w = 360
    var           h = 480

    screenres w, h, 32
   
    var           pixel_size = 4
    
    for y as long = 0 to h step pixel_size
      PixelBits csg, 0, y, pixel_size
    next

    sleep

end sub

Main
Last edited by dafhi on Oct 25, 2018 5:45, edited 1 time in total.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

Post by deltarho[1859] »

@dafhi

I know about LCG's weaknesses but what is your point?
dodicat
Posts: 7979
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

Post by dodicat »

The problem remains, how to get a random ulongint within a given range of ulongints?
Strings can do it, but are slow.
deltarho
your beloved 64 bit (whatever gcc I reckon) is really crap.
So do this with -gen gas first.

Code: Select all



 Function plus(byval num1 As String,byval num2 As String) As String
static As Ubyte AddQMod(0 To 19)={48,49,50,51,52,53,54,55,56,57,48,49,50,51,52,53,54,55,56,57}
static As Ubyte AddBool(0 To 19)={0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1}
        Var n_=0
        Dim As Ubyte addup=Any,addcarry
        #macro finish()
        Return Ltrim(answer,"0")
        #endmacro
        If Len(num2)>Len(num1) Then  Swap num2,num1
        Var diff=Len(num1)-Len(num2)
        Var answer="0"+num1
        For n_=Len(num1)-1 To diff Step -1 
            addup=num2[n_-diff]+num1[n_]-96
            answer[n_+1]=ADDQmod(addup+addcarry)
            addcarry=ADDbool(addup+addcarry)
        Next n_ 
        If addcarry=0 Then 
            finish()
        End If
        If n_=-1 Then 
            answer[0]=addcarry+48
            finish()
            End if
            For n_=n_ To 0 Step -1 
                addup=num1[n_]-48
                answer[n_+1]=ADDQmod(addup+addcarry)
                addcarry=ADDbool(addup+addcarry)
                If addcarry=0 Then Exit For
            Next n_
            answer[0]=addcarry+48
            finish()
        End Function
         Function minus(num1 As String,num2 As String) As String
            Dim  subQmod(0 To 19) As Ubyte
            Dim   subbool(0 To 19) As Ubyte
            For z As Integer=0 To 19
                subqmod(z)=(z Mod 10+48)
                subbool(z)=(-(10>z))
            Next z
            Dim As Integer bigger,swapflag           
            Dim sign As String * 1
            Var lenf=Len(NUM1)
            Var lens=Len(NUM2)
            #macro finishup()
            answer=Ltrim(answer,"0")
            If answer="" Then Return "0"
            If swapflag=1 Then Swap NUM1,NUM2
            Return sign+answer
            #endmacro
            #macro compare()
            If Lens>lenf Then bigger= -1:Goto fin
            If Lens<lenf Then bigger =0:Goto fin
            If NUM2>NUM1 Then 
                bigger=-1
            Else
                bigger= 0
            End If
            fin:
            #endmacro
            
            compare()
            If bigger Then 
                sign="-"
                Swap NUM2,NUM1
                Swap lens,lenf
                swapflag=1
            End If
            Var diff=lenf-lens
            Dim As String answer=NUM1
            Dim As Integer n
            Dim As Ubyte takeaway,subtractcarry
            subtractcarry=0
            For n=lenf-1 To diff Step -1 
                takeaway= num1[n]-num2[n-diff]+10-subtractcarry
                answer[n]=Subqmod(takeaway)
                subtractcarry=Subbool(takeaway)
            Next n 
            
            If subtractcarry=0 Then:finishup():End If
            If n=-1 Then:finishup():End If
            For n=n To 0 Step -1 
                takeaway= num1[n]-38-subtractcarry 
                answer[n]=Subqmod(takeaway)
                subtractcarry=Subbool(takeaway)
                If subtractcarry=0 Then exit for
            Next n
            finishup()
        End Function  
        
  Function Get64Bit() As UlongInt
  Return (Cast( Ulongint, Rnd*(2^32) ) Shl 32) Or Cast( Ulongint, Rnd*(2^32) )   
End Function      
        
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)
    Scope
        Var 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:
    End Scope
    #endmacro
    Dim As Long L=Len(s1),ans=1
    Dim As String s=String(L,0)
    While ans
        GetNumber
        compare(s,s1,ans)
    Wend
    var g=ltrim(s,"0")
    if g="" then g="0"
    Return g
   
End Function

function rndX overload(s1 as ulongint) as ulongint
    static as ulongint u=18446744073709551615
    return iif(s1=u,Get64Bit(),Get64Bit() mod (s1+1))
end function


Type MsWsParams
  Declare Constructor
  Declare Sub MyRandomize( Byval seed As Ulongint = 0, Byval sequence As Ulongint = 0 )
  Declare Function Rand() As Ulongint
  Declare Function RandSE() As Double
  Seed As Ulongint
  Sequence As Ulongint
End Type


Sub MsWsParams.MyRandomize( Byval seed As Ulongint = 0, Byval sequence As Ulongint = 0 )
  Randomize , 5
  If seed = 0 Then
    This.Seed = Get64Bit
    This.Sequence = Get64Bit
  Else
    If sequence = 0 Then
      This.Seed = seed
      This.Sequence = Get64Bit
    Else ' For a sequence repeat
      This.Seed = seed
      This.Sequence = sequence
    End If
  End If
  ' Warm up generator - essential for any PRNG
  For i As Ulong = 1 To 1000
    This.rand()
  Next
End Sub

' 18446744073709551557 Is the largest prime number < 2^64

Function MsWsParams.Rand() As Ulongint
     dim as ulongint x=-1
  This.Seed *= This.Seed : This.Sequence += 18446744073709551557 : This.Seed += This.Sequence
  This.Seed = ( This.Seed Shr 32 ) Or ( This.Seed Shl 32 )
  Return( This.Seed )
End Function

Function MsWsParams.RandSE() As Double
   
  This.Seed *= This.Seed : This.Sequence += 18446744073709551557 : This.Seed += This.Sequence
  
  This.Seed = ( This.Seed Shr 32 ) Or ( This.Seed Shl 32 )
  Return This.Seed/2^64
End Function

Constructor MsWsParams
  This.MyRandomize()
End Constructor

'#Include "MsWs.bas"

Dim shared As MsWsParams MsWs
Dim As Ulong i

'MSWs.MyRandomize(123456, 987654)
'MsWs.MyRandomize( 123456, 0)
dim shared as long outside

randomize ,2 'for string speed
function range overload(byref f as string,byref l as string) as string
    var u= plus (RndX(minus(l,f)),f)
    'for testing range
    if valulng(u)<valulng(f) then u=f:outside+=1
    if valulng(u)>valulng(l) then u=l:outside+=1
    return  u
end function


function range overload(f as ulongint,l as ulongint) as ulongint
dim as ulongint u= Int(MsWs.randse*((l+1)-(f))+(f))
'dim as ulongint u= Int(rnd*((l+1)-(f))+(f))
if u<f then u=f:outside+=1
if u>l then u=l:outside+=1
return u
end function

dim as ulongint k=-1

    outside=0
    k=k\2

 print "ulongints"   
print "Range required"
print k-5;"   to   ";k

redim as long a(0 to 5)

for i=1 to 1000000
    dim as ulongint u=range(k-5,k)
    a(k-u)+=1
next

for n as long=0 to 5
    print a(n)
next
print "done, errors",outside

for n as long=0 to 5
    a(n)=0
    next
print "_____________"
print "strings"

outside=0


print str(k-5);"   to   ";str(k)
for i=1 to 1000000
    dim as string u=range(str(k-5),str(k))
    a((k-valulng(u)) )+=1
next

for n as long=0 to 5
    print a(n)
next
print "done, errors",outside

print "_____________"




sleep
 
 
  
Last edited by dodicat on Oct 25, 2018 2:24, edited 2 times in total.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

Post by deltarho[1859] »

dodicat wrote:your beloved 64 bit (whatever gcc I reckon) is really crap.
Can you be more specific? What beloved 64 bit?
dodicat
Posts: 7979
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

Post by dodicat »

fbc 1.06.0 64 bit with gcc 8.1/8.2 (maybe)
Or 1.05.0 64 bit fbc.(possibly)
I thought I read recently that you are sticking with 64 bits (somebody else , forgot who, said that they wouldn't stoop to even run a test on the 32 bit compiler)
Such is life I suppose.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

Post by deltarho[1859] »

From the opening post:
All of those figures are based upon 9 runs each and the median taken for FBC 1.05/gcc 5.2 64-bit.
BTW, I don't do beloved anything.
dodicat
Posts: 7979
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

Post by dodicat »

It's alright deltarho.
You are a reasonably level headed guy ( PractRand gold medal split for weight distribution).
A chip on each shoulder does the same job for me.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

Post by deltarho[1859] »

I have been using this for ages - why hasn't anyone pulled me on it?

Code: Select all

Function Get64Bit() As Ulongint
  Return (Cast( Ulongint, Rnd*(2^32) ) Shl 32) Or Cast( Ulongint, Rnd*(2^32) )
End Function

Code: Select all

#Define Get64Bit Cast( Ulongint, Rnd*2^64 )
dafhi
Posts: 1641
Joined: Jun 04, 2005 9:51

Re: Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

Post by dafhi »

@deltarho - the point is that luck isn't a factor.

anyway, again, thanks for this cool framework
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

Post by deltarho[1859] »

dafhi wrote:@deltarho - the point is that luck isn't a factor.
I am sorry but I have just lost interest in talking with you.
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

Post by jj2007 »

dodicat wrote:The problem remains, how to get a random ulongint within a given range of ulongints?
Like this, maybe?

Code: Select all

Dim As ulongint MyMin=100, MyMax=200, MyRand
Dim As Double MyRange=(MyMax-MyMin)
  for i As integer=1 to 500
	MyRand=Rnd()*MyRange+MyMin
	Print MyRand; " ";
  next
  sleep
Post Reply