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

General FreeBASIC programming questions.
deltarho
Posts: 2003
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 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 on Nov 03, 2018 9:30, edited 12 times in total.
dafhi
Posts: 1251
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: 2003
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 = 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: 1242
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: 1251
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 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
Posts: 2003
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: 5947
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 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
#macro finish()
#endmacro
If Len(num2)>Len(num1) Then  Swap num2,num1
Var diff=Len(num1)-Len(num2)
For n_=Len(num1)-1 To diff Step -1
Next n_
finish()
End If
If n_=-1 Then
finish()
End if
For n_=n_ To 0 Step -1
If addcarry=0 Then Exit For
Next n_
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()
If answer="" Then Return "0"
If swapflag=1 Then Swap NUM1,NUM2
#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
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
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 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
Posts: 2003
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: 5947
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: 2003
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: 5947
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: 2003
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: 1251
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: 2003
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: 1242
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, MyRand
Dim As Double MyRange=(MyMax-MyMin)
for i As integer=1 to 500
MyRand=Rnd()*MyRange+MyMin
Print MyRand; " ";
next
sleep