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

General FreeBASIC programming questions.
deltarho
Posts: 1718
Joined: Jan 02, 2017 0:34
Location: UK

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

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 UlongintEnd TypeSub 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()  NextEnd Sub' 18446744073709551557 is the largest prime number < 2^64Function 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 FunctionFunction 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^64End FunctionFunction 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 dafhiFunction MsWsParams.Range(First As Longint, Last As Longint) As LongInt  Function = This.Seed/2^64 * (Last - First + 1) - .5 + FirstEnd FunctionConstructor MsWsParams  This.MyRandomize()End Constructor`
Last edited by deltarho on Nov 03, 2018 9:30, edited 12 times in total.
dafhi
Posts: 1234
Joined: Jun 04, 2005 9:51

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

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
Posts: 1718
Joined: Jan 02, 2017 0:34
Location: UK

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

@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 = statestate = mul * state + add ' honours generatoroutput = 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: 1092
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

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

deltarho 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: 1234
Joined: Jun 04, 2005 9:51

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

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 functionconst             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    NextEnd Subsub 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    sleepend subMain`
Last edited by dafhi on Oct 25, 2018 5:45, edited 1 time in total.
deltarho
Posts: 1718
Joined: Jan 02, 2017 0:34
Location: UK

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

@dafhi

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

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

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 Stringstatic 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=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=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=range(48,s1)    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 Functionfunction rndX overload(s1 as ulongint) as ulongint    static as ulongint u=18446744073709551615    return iif(s1=u,Get64Bit(),Get64Bit() mod (s1+1))end functionType 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 UlongintEnd TypeSub 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()  NextEnd Sub' 18446744073709551557 Is the largest prime number < 2^64Function 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 FunctionFunction 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^64End FunctionConstructor MsWsParams  This.MyRandomize()End Constructor'#Include "MsWs.bas"Dim shared As MsWsParams MsWsDim As Ulong i'MSWs.MyRandomize(123456, 987654)'MsWs.MyRandomize( 123456, 0)dim shared as long outsiderandomize ,2 'for string speedfunction 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  uend functionfunction range overload(f as ulongint,l as ulongint) as ulongintdim 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+=1if u>l then u=l:outside+=1return uend functiondim as ulongint k=-1    outside=0    k=k\2 print "ulongints"   print "Range required"print k-5;"   to   ";kredim as long a(0 to 5)for i=1 to 1000000    dim as ulongint u=range(k-5,k)    a(k-u)+=1nextfor n as long=0 to 5    print a(n)nextprint "done, errors",outsidefor n as long=0 to 5    a(n)=0    nextprint "_____________"print "strings"outside=0print 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)) )+=1nextfor n as long=0 to 5    print a(n)nextprint "done, errors",outsideprint "_____________"sleep    `
Last edited by dodicat on Oct 25, 2018 2:24, edited 2 times in total.
deltarho
Posts: 1718
Joined: Jan 02, 2017 0:34
Location: UK

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

dodicat wrote:your beloved 64 bit (whatever gcc I reckon) is really crap.

Can you be more specific? What beloved 64 bit?
dodicat
Posts: 5614
Joined: Jan 10, 2006 20:30
Location: Scotland

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

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
Posts: 1718
Joined: Jan 02, 2017 0:34
Location: UK

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

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: 5614
Joined: Jan 10, 2006 20:30
Location: Scotland

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

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
Posts: 1718
Joined: Jan 02, 2017 0:34
Location: UK

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

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: 1234
Joined: Jun 04, 2005 9:51

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

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

anyway, again, thanks for this cool framework
deltarho
Posts: 1718
Joined: Jan 02, 2017 0:34
Location: UK

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

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: 1092
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

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

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, MyRandDim As Double MyRange=(MyMax-MyMin)  for i As integer=1 to 500   MyRand=Rnd()*MyRange+MyMin   Print MyRand; " ";  next  sleep`