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

General FreeBASIC programming questions.
MrSwiss
Posts: 3222
Joined: Jun 02, 2013 9:27
Location: Switzerland

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

dodicat wrote:It is also much faster than dipping into double and multiplying.
The uniformity seems OK. Perhaps the randomness needs further testing.

Agreed, but only (IMO), if you can use the generated value unmodified.
Putting it into a range, creates the problems ...
deltarho[1859]
Posts: 1930
Joined: Jan 02, 2017 0:34
Location: UK

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

Here is a very easy way to test the validity of an integer range function instead of dividing the range into equal units and counting how many 'hits' the units get looking for a reasonable balance. If we have a uniform distribution then the average will be equal to (first + last)/2 for an infinite test but fairly close with a 'decent' size test.

Code: Select all

`Function rnd_range (first As Longint, last As Longint) As LongInt  Function = Rnd*(last - first + 1) - .5 + firstEnd Function Randomize Dim As Ulong iDim As Longint tot For i = 1 To 10^8  tot += rnd_range(-123,456)NextPrint tot/10^8 Sleep`

Examples:

Code: Select all

`(-5,5)     =>   0    Progamatically I got 1.118e-005(-10,-1)   =>  -5.5   -do- I got -5.49982073(0,9)      =>   4.5   -do- I got 4.49988167(1,10)     =>   5.5   -do- I got 5.50016871(347,692)  => 519.5   -do- I got 519.50917841(-123,456) => 166.5   -do- I got 166.48569356`

That very quickly proved that '-(first>=0)' was actually wrong.

@dodicat

Still not entirely sure what you mean but will this do? It is simply the above adapted to take a random 64-bit value as a further parameter.

Code: Select all

`Function rnd_range (RandomNumber as UlongInt, first As Longint, last As Longint) As LongInt  Function = RandomNumber/2^64*(last - first + 1) - .5 + firstEnd Function`
dodicat
Posts: 5913
Joined: Jan 10, 2006 20:30
Location: Scotland

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

similar results with no float operators, only uongint/longint operations.

Code: Select all

` '==========================type Rand   a as ulongint   b as ulongint   c as ulongint   d as ulongintend typefunction ranulong(x as rand) as ulongint   dim e as ulongint = x.a - ((x.b shl 7) or (x.b shr (64 - 7)))   x.a = x.b xor ((x.c shl 13) or (x.c shr (64 - 13)))   x.b = x.c + ((x.d shl 37) or (x.d shr (64 - 37)))   x.c = x.d + e   x.d = e + x.a   return x.dend function'function randouble(x as rand) as double ''NOT NEEDED    'return ranulong(x)/18446744073709551615ull   ' end function sub init(x as rand, byval seed as ulongint=100)   dim i as ulongint    x=type(4058668781,seed,seed,seed)    for i as ulongint=1 to 20        ranulong(x)        nextend sub'=========dim shared as rand zinit(z,200)'========function range(f as longint,l as longint) as longint 'for any ranges, ulongint for + only ranges    return (ranulong(z) mod ((l-f)+1)) + fend functionDim As Ulong iDim As Longint tot dim as double t=timer For i = 1 To 10^8  tot += range(-5,5)NextPrint tot/10^8tot=0For i = 1 To 10^8  tot += range(-10,-1)NextPrint tot/10^8tot=0For i = 1 To 10^8  tot += range(0,9)NextPrint tot/10^8tot=0For i = 1 To 10^8  tot += range(1,10)NextPrint tot/10^8tot=0For i = 1 To 10^8  tot += range(347,692)NextPrint tot/10^8tot=0For i = 1 To 10^8  tot += range(-123,456)NextPrint tot/10^8print "Time taken  ";timer-tsleep     `

Code: Select all

` -0.00051191-5.50043617 4.49987626 5.50025758 519.49842485 166.4962325Time taken   2.242408916354179 `
sean_vn
Posts: 283
Joined: Aug 06, 2012 8:26

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

Nice RNG algorithm. I must remember it.
There is a 128 bit multiply in the AMD 64 instruction set that you can use for fixed point calculations. A random 64 bit number by (n+1) will give a result between 0 and n in the register that contains the high 64 bits of the 128 bit result. Not perfect but generally okay.
deltarho[1859]
Posts: 1930
Joined: Jan 02, 2017 0:34
Location: UK

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

@dodicat

You are using a few rotations there and they remind me of Sebastiano Vigna's xorshiro algorithms where he used a C macro

Code: Select all

`static inline uint64_t rotl(const uint64_t x, int k) {  return (x << k) | (x >> (64 - k));}`

I still haven't got my head around

Code: Select all

`Return (ranulong(z) Mod ((l-f)+1)) + f`

yet but I will get there. <smile>

It does not look like it will compete against Knuth64 on speed but I bet it has better randomness than Knuth64 which leaves a lot to be desired.
dodicat
Posts: 5913
Joined: Jan 10, 2006 20:30
Location: Scotland

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

deltarho
The difference is ( last-first +1)
then first + random number mod difference puts the final number in the output range.
The randomising algo can be half decent ulongint generator.
The range function simply avoids doubles.
The mod operator works in the ulongint range, although the help file says integer.
You can overload for longint and ulongint.
pass parameters with ll or ull subscripts.
The line
return (ranulong(z) mod ((l-f)+1)) + f
is generally
return (someulongintgenerator(z) mod ((l-f)+1)) + f
You own generator would be very fast I think.
deltarho[1859]
Posts: 1930
Joined: Jan 02, 2017 0:34
Location: UK

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

dodicat wrote:You own generator would be very fast I think.

Blindingly fast but it failed the fast validity test. [1]

I used this in the opening post:

Code: Select all

`Function MsWsParams.Range( First As Longint, Last As Longint ) As Longint  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 Mod ( Last - First + 1)  + FirstEnd Function`

Added: [1]. For the life of me, I cannot see how Mod can possibly work in this context.
dodicat
Posts: 5913
Joined: Jan 10, 2006 20:30
Location: Scotland

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

For your generator I use

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    Declare Function Range Overload ( First As Ulongint, Last As Ulongint ) As Ulongint    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'alterationsFunction MsWsParams.Range( First As Longint, Last As Longint ) As Longint    Return (this.Rand() Mod ((last-first)+1)) + first End FunctionFunction MsWsParams.Range( First As Ulongint, Last As Ulongint ) As Ulongint    Return (this.Rand() Mod ((last-first)+1)) + first End FunctionConstructor MsWsParamsThis.MyRandomize()End ConstructorDim  As MsWsParams z Dim As Ulong iDim As Longint totDim As Double t=Timer For i = 1 To 10^8    tot += z.range(-5ll,5)NextPrint tot/10^8tot=0For i = 1 To 10^8    tot += z.range(-10ll,-1)NextPrint tot/10^8tot=0For i = 1 To 10^8    tot += z.range(0ull,9)NextPrint tot/10^8tot=0For i = 1 To 10^8    tot += z.range(1ull,10)NextPrint tot/10^8tot=0For i = 1 To 10^8    tot += z.range(347ull,692)NextPrint tot/10^8tot=0For i = 1 To 10^8    tot += z.range(-123ll,456)NextPrint tot/10^8Print "pure ulongint example"For i = 1 To 10    Print z.range(18446744073709551610ull,18446744073709551613)NextPrint "Time taken  ";Timer-tSleep  `

result

Code: Select all

`-6.353e-005-5.50021221 4.50014817 5.49993317 519.49835457 166.50732974pure ulongint example18446744073709551611184467440737095516101844674407370955161318446744073709551613184467440737095516101844674407370955161218446744073709551612184467440737095516111844674407370955161218446744073709551613Time taken   2.048540160991252 `
deltarho[1859]
Posts: 1930
Joined: Jan 02, 2017 0:34
Location: UK

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

@dodicat

I was using This.Seed/2^64 when I should have used just This.Seed. Dear, oh dear. Your method is producing the same results as dafhi's and my Int method before that.

The reason why I used This.Seed, as opposed to This.Rand, is because the engine only takes up a few lines so I use the engine for every function to avoid calling another function. As it turns out the difference between the engine plus This.Seed and This.Rand is negligible. FB keeps surprising me.

Anyway, my original Int method, using '-gen gcc -Wc -O3' in 32-bit mode, was beaten by dafhi's method from 34MHz to 63MHz. Your method storms in at 135MHz.

The most extraordinary thing is using '-gen gcc -Wc -O3' in 64-bit mode saw my Int method and dafhi's method coming in at the same speed and now your method is the same. Goodness knows what the binary looks like and I find it hard to believe that all three methods have been reduced to the same native code. However, they are coming in at the same speed nonetheless.

Since you MODification has improved the 32-bit results I have revised the LongInt in the opening post yet again for folks who use 32-bit.

I must now sit down and fathom out how the Mod actually works. How did you come up with that? Is it from your archives or did you come up with it during the life of this thread? This is important because if you got it from your archives then you qualify for a fortnight in Clacton-On-sea. On the other hand, if you came up with it during the life of this thread you qualify for a week in Clacton-On-Sea. The burning question is are you going to tell 'porkies' to avoid going to Clacton for a fortnight.

I wonder what MrSwiss thinks about your Mod approach. We cannot get much more elegant than that and it is very close to being four times faster than my original Int method. That should definitely go into the manual but it will, I am sure, leave a lot of folks scratching their heads.

Actually, I am being unfair about Clacton because I have stayed in a few Bed & Breakfasts there and I have no complaints on that score. It is probably very good for families as there is plenty of stuff going on to keep the kids fully occupied. It is over 25 years since I was last there but there was a superb fish restaurant that I used to frequent owned by a divorced lady who I took a shine too. She knocked out a haddock to die for. <smile>
deltarho[1859]
Posts: 1930
Joined: Jan 02, 2017 0:34
Location: UK

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

Remembering that this thread is about MsWs it is worth noting that in 64-bit mode MsWs' range function is about 11% faster than PCG32II. However, for RandSE PCG32II is about 12% faster than MsWs. In 'real life' scenarios it may not be easy to separate them. There is a difference in that PCG32II is a 32-bit generator and MsWs is a 64-bit generator with PCG32II's RandSE having 32-bit granularity and MsWs' RandSE having 53-bit granularity. It could be argued then that MsWs has elbowed PCG32II off the podium.
dodicat
Posts: 5913
Joined: Jan 10, 2006 20:30
Location: Scotland

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

My ideal fortnight would be absolute silence and solitude where I can practice my medication techniques, thus I would choose Morecambe.
MrSwiss
Posts: 3222
Joined: Jun 02, 2013 9:27
Location: Switzerland

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

dodicat, just make very shure to NOT run out of medicine or have a established supply-chain there.<lol>
dafhi
Posts: 1245
Joined: Jun 04, 2005 9:51

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

range using modulus may fit unevenly into 2^64. the remainder goes to the least significant digits
deltarho[1859]
Posts: 1930
Joined: Jan 02, 2017 0:34
Location: UK

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

Code: Select all

`Dim As Longint lo = 20, hi = 147, R, ansDim As Double xR = 7547257594037927936ans =  R Mod (hi - lo + 1)  + 20 ' dodicatPrint ansx = R/2^64*(hi - lo + 1) - 0.5 + 20ans = R/2^64*(hi - lo + 1) - 0.5 + 20  ' dafhiPrint ans, xSleep`

gives

Code: Select all

` 20 72            71.86961971048731`

dafhi's ans is correct

Perhaps we need to write our own Mod function if it is fast enough.
dafhi
Posts: 1245
Joined: Jun 04, 2005 9:51

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

also, you could get rid of the hi (+ 1) if you change params to double

deltarho[1859] wrote:Perhaps we need to write our own Mod function if it is fast enough.

:-)

i'm having a blast with this RNG stuff

Return to “General”

### Who is online

Users browsing this forum: No registered users and 5 guests