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

[dodicat]

Thanks jj2007.
I have used your code, but adjusted the myrand with int(~~) , to give a more uniform distribution, which is shown.
Also, I have noticed a little error in my previous range functions, (the positions of the brackets within)

My gripe is not with the mode of getting rnd(), but with the way gcc seems to fail within high ranges.
This seems OK with -gen gas, but fails with -gen gcc, try with 64 bit fbc for example.

Code: Select all

  
  
  dim as ulongint k=-1
  k=k\2-5
                      '4294967295            ulong max
                      '18446744073709551615  ulongint max
dim as ulongint offset=12345678909876543210
Dim As ulongint MyMin=0, MyMax=5, MyRand

mymin=offset+mymin
mymax=offset+mymax

print mymin;"  to  ";mymax
dim as ulong a(0 to mymax - mymin)

dim as ulongint errors
Dim As ulongint MyRange=(MyMax-MyMin)
  for i As integer=1 to 10000000
   ''MyRand=Rnd()*MyRange+MyMin
    MyRand=int(Rnd()*(MyRange+1))+MyMin
    if myrand>mymax then errors+=1
    if myrand<mymin then errors+=1
   a(mymax-myrand)+=1
  next
  
  print "distribution"
  for n as long=lbound(a) to ubound(a)
      print n, a(n)
  next
  print "errors ";errors
  sleep
   

Of course I could be not seeing something obvious.

[jj2007]

Seems to be a problem with Rnd()'s internal conversion from int to double; try if myrand<mymin then errors-=1, you will see that this never gets triggered. The upper end is the problem.

Remember that most random generators natively return an integer, not a double. If you really need an integer, go to the source, e.g. of PCG32 ;-)

[dodicat]

I can set mymin to 10 and mymax to 15 say.
So I have to check each end.

Anyway, I'll have a mess around from a source (ulongint), and stop using rnd*

[MrSwiss]

Quick & dirty ULongInt generator, using Mersenne Twister ...

Code: Select all

' may NOT work (without warnings) in FBC 32
' only FBC 64 tested!
Randomize(Timer, 3) ' Mersenne Twister (FB internal)
For n As UInteger = 1 To 50 ' 300 call's (warm up)
    Rnd() : Rnd() : Rnd() : _
    Rnd() : Rnd() : Rnd()
Next

Function rand64() As ULongInt           ' uses MT
    Return ((Rnd() * &hFFFFFFFFul) Shl 32) + (Rnd() * &hFFFFFFFFul)
End Function

' quick test
For i As UInteger = 1 To 20
    ? Tab(3); rand64()
Next

? : ? : ? "press a key to EXIT ... ";
Sleep
' ----- EOF -----

[dodicat]

Gcc seems to be a problem here.
Here is a generator without rnd().
It passed deltarho's random test (practrand).

Code: Select all

 


type Rand
	a as ulongint
	b as ulongint
	c as ulongint
	d as ulongint
end type

function 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.d
end function

function randouble(x as rand) as double
    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)
        next
end sub

'=========
dim as rand z
init(z)
'========

 dim as ulongint k=-1
  k=k\2-5
                      '4294967295            ulong max
                      '18446744073709551615  ulongint max
dim as ulongint offset=12345678909876543210
Dim As ulongint MyMin=10, MyMax=15, MyRand

mymin=offset+mymin
mymax=offset+mymax

print mymin;"  to  ";mymax
dim as ulong a(0 to mymax - mymin)

dim as ulongint errors
Dim As ulongint MyRange=(MyMax-MyMin)
  for i As integer=1 to 1000000
   ''MyRand=Rnd()*MyRange+MyMin
    MyRand=int(randouble(z)*(MyRange+1))+MyMin
    if myrand>mymax then myrand=mymax:errors+=1
    if myrand<mymin then myrand=mymin:errors+=1
   a(mymax-myrand)+=1
  next
  
  print "distribution"
  for n as long=lbound(a) to ubound(a)
      print n, a(n)
  next
  print "errors ";errors
  sleep
 

 

[dodicat]

I tested a ulonint range with 64 bit pascal (qword~ulongint).
similar results to -gen gcc freebasic.

Code: Select all

64 bit pascal
8234567890987664  to  8234567890987669
distribution
0    1666458
1    1668248
2    1666887
3    1665010
4    1666652
5    1666745
errors  0



64 bit freebasic
8234567890987664  to  8234567890987669
distribution
 0            1664939
 1            1667938
 2            1666083
 3            1667327
 4            1666684
 5            1667029
errors 0

So -gen gas seems to be a bonus in this context.(can use all of ulongint range)
So it is up to the particular longint algo for randomness within the range.

Deltarho;
I had made a bloomer with your random.
the range function should be(I reckon)

Code: Select all

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

This gives a good distribution and no errors.
Sorry.

[deltarho[1859]]

Gentlemen, it is not the Ulongint range at issue but significant figure accuracy. Remember Ulongint is 19+ and Double is 15+. Push those limits and we should expect trouble.

Many years ago Jim Wilkinson FRS, one of my Numerical Analysis heroes, was at the National Physical Laboratory (NPL) here in the UK and the powers that be wanted to scrap Colossus, a behemoth of a computer and replace it with up to date kit. He threatened to resign on the grounds that the up to date kit had nothing like the significant figure accuracy of the ancient Colossus. He won and Colossus survived a few more years. I met Wilkinson at Colchester University over 40 years ago when he was a guest speaker at a conference my Ministry of Defence sent me to. Wilkinson died in 1986 at Teddington where NPL is based.

[deltarho[1859]]

deleted

[deltarho[1859]]

I have added a couple of overloaded range functions to the opening post: One is Double based and the other is Ulongint based.

The opening post refers to 2nd revision:

Original update used 'Function = This.Seed/2^64*( Last - First ) + First' for the Ulongint version. For a range of 1 to 10, for example, this gave a poor distribution for 1 and 10.
Revised to 'Function = Int(This.Seed/2^64*( Last - First + 1) + First)' gives a balanced distribution.

Usage:

Code: Select all

#Include "MsWs.bas"
Dim as MsWsParams MsWs
 
Dim as Ulong i
For i = 1 To 16
  Print MsWs.Range( 1, 10 )
Next
Print
For i = 1 to 16
  Print MsWs.Range( 1., 2. )
Next
Sleep

Typical output:

Code: Select all

5
7
7
10
8
10
6
5
7
2
4
2
1
8
10
9

 1.911959986558758
 1.372155961150833
 1.720977375183379
 1.186729199672825
 1.95335062271225
 1.086848098707967
 1.134831185547908
 1.256109807795861
 1.158826337460855
 1.274496901540458
 1.826109401664264
 1.155801271430077
 1.949176785107129
 1.413342691679567
 1.488260385432489
 1.754680250121552

[MrSwiss]

Revised to 'Function = Int(This.Seed/2^64*( Last - First + 1) + First)' gives a balanced distribution.

The disribution problem relates to: Int()'s rounding preference.
This is the reason that I'm only using: CInt() (never Int()) but this creates another problem:
Integer = a Long in FBC 32, or a LongInt in FBC 64.

You've written that it should be: ULongInt (or LongInt?) anyway, a different conversion is needed:
'Function = CULngInt(This.Seed/2^64*( Last - First) + First)' ... (or CLngInt()?)

[deltarho[1859]]

The disribution problem relates to: Int()'s rounding preference.

No, the problem relates to a Function returning UlongInt.

If a value is determined as 12.4, for example, then 12 will be returned by the function. On the other hand, if a value is determined as 12.6 then 13 will be returned by the function.

So, if we consider a range of 1 to 10 and store their count, requesting a large number, in a(1 to 10) half of the values, on average, less than 1 will increment a(2). Between 1 and 2 will see half of them, on average, increment a(3). For values which are greater or equal to 10 will see half them increment a non-existent a(11).

For 100 million requests we get a distribution like this:

Code: Select all

5554481
11112826
11110383
11109798
11113737
11112768
11111364
11111061
11108670
5554912

a(2) to a(9) look OK but they have built incorrectly. a(7), for example, is getting, on average, half its count from a(6) and passing half of its own count, on average, to a(8).

By adding 1 to Last - First and using Int ( Returns the floor of a number ) we get a distribution like this:

Code: Select all

9998252
10002953
10001378
9998464
9997937
10003603
10000452
9999618
9996366
10000977

'Function = CULngInt(This.Seed/2^64*( Last - First) + First)' will also round up, given half the chance [half the chance - get it? Never mind] and we will get a distribution along the lines of the first one above.

So, 'Function = Int( This.Seed/2^64*( Last - First + 1 ) + First )' is one way to get a balanced distribution. Whether we compile to 32-bit or 64-bit is neither here nor there. We do not get an issue with the overloaded Double range.

[deltarho[1859]]

I don't think that I have ever used an unsigned range and why I used UlongInt. However, some folk might so I have changed UlongInt to LongInt.

The Function return has been changed to:

Code: Select all

Function = Int( This.Seed/2^64*( Last - First - (First>=0) ) + First )

If First>=0 then it follows that Last>0 so we use 'Last - First + 1' as we did when using UlongInt.

However, if First<0 then we use 'Last - First' and do not add 1.

Opening post now refers to 3rd revision.

Bit of a can of worms - I didn't expect this. <smile>

[MrSwiss]

Bit of a can of worms - I didn't expect this. <smile>

Seems to be a matter of underestimating the problem.

Btw: Int() produces at maximum, a Double with 31 bit granularity (bad news!). FBC 32.
No such problem with FBC 64.

The better way, to deal with it is IMO, to resolve rounding of ULongInt (instead of using Int()).
As done below:

Code: Select all

Randomize(Timer, 3)

For n As UInteger = 1 To 100
    Rnd() : Rnd() : Rnd()
Next
'===========================

Dim As ULongInt a(0 To 10), t, ll = 0, hl = 10, lcnt, hcnt, cnt

For z As UInteger = 1 To 100
    For i As ULongInt = 0 To 1e6 - 1
        again:
        cnt += 1
        t = CULngInt( Rnd * ((hl + 0.5) - (ll - 0.5)) + (ll - 0.5) )
        If t < ll Then
            lcnt += 1 : GoTo again
        ElseIf t > hl Then
            hcnt += 1 : GoTo again
        Else
            a(t) += 1
        End If
    Next
next

For i As UInteger = 0 To 10
    Print i; Tab(9); a(i)
Next

Print
Print "to high: "; hcnt; "   to low: "; lcnt; "   count: "; cnt
Print "out of range: "; (hcnt + lcnt)
Print
Print "press a key to EXIT ... "
Sleep

[deltarho[1859]]

Seems to be a matter of underestimating the problem.

You mean as you did when you suggested CULngInt. <Big Grin>

Btw: Int() produces at maximum, a Double with 31 bit granularity (bad news!). FBC 32.

Can you give an example - I cannot prove that statement.

As done below:

On an elegance scale of 10: MrSwiss 1; Deltarho 7

Can anyone come up with an 8 or more?

Prizes will be given for the best entries: 1st Prize - A week in Clacton-On-Sea; 2nd Prize - A fortnight in Clacton-on-Sea. Don't get it? Visit Clacton-On-Sea and the penny will drop.

[MrSwiss]

On an elegance scale of 10: MrSwiss 1; Deltarho 7

On an efficency scale of 10: MrSwiss 10; Deltarho 3 (none ?? duno)

I don't code, to win elegance competitions (I want correctly working code).

You mean as you did when you suggested CULngInt.

What's your case? That problem was easily solved! (with inelegant method)