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:

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

Post by deltarho[1859] »

dafhi wrote:also, you could get rid of the hi (+ 1) if you change params to double
and the -0.5

Double has been taken care of.

In any case, we should focus on one thing at a time. Mod is the issue.
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:The mod operator works in the ulongint range, although the help file says integer.
'R Mod n' returns zero if R >= 2^32-1.

If we write our own Mod function it looks like we need extended precision, ie 80-bit, but we don't have it. We should have since the FPU uses 80-bit. PowerBASIC has extended precision.
MrSwiss
Posts: 3910
Joined: Jun 02, 2013 9:27
Location: Switzerland

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

Post by MrSwiss »

deltarho[1859] wrote:'R Mod n' returns zero if R >= 2^32-1.
This sort of statement is the reason, that I don't trust code of yours.

Proof that above statement is nonsense (tested FBC 32/64, 1.05.0):

Code: Select all

Const As LongInt    pllmax = &h7FFFFFFFFFFFFFFFll, _    ' positive 64 bit maximum
                    nllmax = &hFFFFFFFFFFFFFFFFll       ' negative 64 bit maximum


Print nllmax Mod 19
Print pllmax Mod 19

Print : Print "done ... ";

Sleep
result wrote:-1
17
Added: the neg. max. referres to binary = set bit's (value -1)
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

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

Post by dodicat »

mod seems OK for > 2^32-1
Why do you reckon it is not?

Here are the string mods and fb mods for ulongints high in the range:

Code: Select all

 
 


Function _divide(n1 As String,n2 As String,decimal_places As integer=10,dpflag As String="s") As String
          Dim As String number=n1,divisor=n2
          dpflag=lcase(dpflag)
          'For MOD
          dim as integer modstop
          if dpflag="mod" then 
              if len(n1)<len(n2) then return n1
              if len(n1)=len(n2) then
                  if n1<n2 then return n1
                  end if
              modstop=len(n1)-len(n2)+1
              end if
          if dpflag<>"mod" then
     If dpflag<>"s"  Then dpflag="raw" 
     end if
        Dim runcount As integer
        '_______  LOOK UP TABLES ______________
        Dim Qmod(0 To 19) As Ubyte
        Dim bool(0 To 19) As Ubyte
        For z As Integer=0 To 19
    Qmod(z)=(z Mod 10+48)
    bool(z)=(-(10>z))
Next z
Dim answer As String   'THE ANSWER STRING  

'_______ SET THE DECIMAL WHERE IT SHOULD BE AT _______
Dim As String part1,part2
#macro set(decimal)
#macro insert(s,char,position)
If position > 0 And position <=Len(s) Then
part1=Mid$(s,1,position-1)
part2=Mid$(s,position)
s=part1+char+part2
End if
#endmacro
insert(answer,".",decpos)
  answer=thepoint+zeros+answer
If dpflag="raw" Then
    answer=Mid(answer,1,decimal_places)
    End if
#endmacro
'______________________________________________
'__________ SPLIT A STRING ABOUT A CHARACTRR __________
Dim As String var1,var2
    Dim pst As integer
      #macro split(stri,char,var1,var2)
    pst=Instr(stri,char)
    var1="":var2=""
    If pst<>0 Then
    var1=Rtrim(Mid(stri,1,pst),".")
    var2=Ltrim(Mid(stri,pst),".")
Else
    var1=stri
    End if
    #endmacro
    
       #macro Removepoint(s)
       split(s,".",var1,var2)
#endmacro
'__________ GET THE SIGN AND CLEAR THE -ve __________________
Dim sign As String
          If Left(number,1)="-" Xor Left (divisor,1)="-" Then sign="-"
            If Left(number,1)="-" Then  number=Ltrim(number,"-")
            If Left (divisor,1)="-" Then divisor=Ltrim(divisor,"-")
              
'DETERMINE THE DECIMAL POSITION BEFORE THE DIVISION
Dim As integer lennint,lenddec,lend,lenn,difflen
split(number,".",var1,var2)
lennint=Len(var1)
split(divisor,".",var1,var2)
lenddec=Len(var2)

If Instr(number,".") Then 
    Removepoint(number)
    number=var1+var2
    End if
If Instr(divisor,".") Then 
    Removepoint(divisor)
    divisor=var1+var2
    End if
Dim As integer numzeros
numzeros=Len(number)
number=Ltrim(number,"0"):divisor=Ltrim (divisor,"0")
numzeros=numzeros-Len(number)
lend=Len(divisor):lenn=Len(number)
If lend>lenn Then difflen=lend-lenn
Dim decpos As integer=lenddec+lennint-lend+2-numzeros 'THE POSITION INDICATOR
Dim _sgn As Byte=-Sgn(decpos)
If _sgn=0 Then _sgn=1
Dim As String thepoint=String(_sgn,".") 'DECIMAL AT START (IF)
Dim As String zeros=String(-decpos+1,"0")'ZEROS AT START (IF) e.g. .0009
if dpflag<>"mod" then
If Len(zeros) =0 Then dpflag="s"
end if
Dim As integer runlength
If Len(zeros) Then 
     runlength=decimal_places
     answer=String(Len(zeros)+runlength+10,"0")
    If dpflag="raw" Then 
        runlength=1
        answer=String(Len(zeros)+runlength+10,"0")
        If decimal_places>Len(zeros) Then
            runlength=runlength+(decimal_places-Len(zeros))
            answer=String(Len(zeros)+runlength+10,"0")
            End If
            End If

Else
decimal_places=decimal_places+decpos
runlength=decimal_places
answer=String(Len(zeros)+runlength+10,"0")
End if
'___________DECIMAL POSITION DETERMINED  _____________

'SET UP THE VARIABLES AND START UP CONDITIONS
number=number+String(difflen+decimal_places,"0")
        Dim count As integer
        Dim temp As String
        Dim copytemp As String
        Dim topstring As String
        Dim copytopstring As String
        Dim As integer lenf,lens
        Dim As Ubyte takeaway,subtractcarry
        Dim As integer n3,diff
       If Ltrim(divisor,"0")="" Then Return "Error :division by zero"   
        lens=Len(divisor)
         topstring=Left(number,lend)
         copytopstring=topstring
        Do
            count=0
        Do
            count=count+1
            copytemp=temp
    
            Do
'___________________ QUICK SUBTRACTION loop _________________              
            
lenf=Len(topstring)
If  lens<lenf=0 Then 'not
If Lens>lenf Then
temp= "done"
Exit Do
End if
If divisor>topstring Then 
temp= "done"
Exit Do
End if
End if

  diff=lenf-lens
        temp=topstring
        subtractcarry=0
        
        For n3=lenf-1 To diff Step -1
            takeaway= topstring[n3]-divisor[n3-diff]+10-subtractcarry
            temp[n3]=Qmod(takeaway)
            subtractcarry=bool(takeaway)
        Next n3 
        If subtractcarry=0 Then Exit Do
         If n3=-1 Then Exit Do
        For n3=n3 To 0 Step -1 
            takeaway= topstring[n3]-38-subtractcarry
             temp[n3]=Qmod(takeaway)
            subtractcarry=bool(takeaway)
            if subtractcarry=0 then exit do
            Next n3
        Exit Do
        
        Loop 'single run
        temp=Ltrim(temp,"0")
        If temp="" Then temp= "0"
            topstring=temp
        Loop Until temp="done"
     ' INDIVIDUAL CHARACTERS CARVED OFF ________________       
        runcount=runcount+1
       If count=1 Then
           topstring=copytopstring+Mid(number,lend+runcount,1)
           Else
       topstring=copytemp+Mid(number,lend+runcount,1)
   End If
       copytopstring=topstring
       topstring=Ltrim(topstring,"0")
       if dpflag="mod" then
       if runcount=modstop then 
           if topstring="" then return "0"
           return mid(topstring,1,len(topstring)-1)
           end if
       end if
       answer[runcount-1]=count+47
       If topstring="" And runcount>Len(n1)+1 Then
           Exit Do
           End if
   Loop Until runcount=runlength+1
   
   ' END OF RUN TO REQUIRED DECIMAL PLACES
   set(decimal) 'PUT IN THE DECIMAL POINT
  'THERE IS ALWAYS A DECIMAL POINT SOMEWHERE IN THE ANSWER
  'NOW GET RID OF IT IF IT IS REDUNDANT
       answer=Rtrim(answer,"0")
       answer=Rtrim(answer,".")
       answer=Ltrim(answer,"0")
       If answer="" Then Return "0"
   Return sign+answer
End Function

#define mod_(a,b) _divide(a,b,,"mod")
#define div_(a,b) _divide(a,b,-2)

dim as ulongint R=2^32-1   + 1000000
dim as ulongint n=234
print R Mod n' returns zero if R >= 2^32-1.

dim as ulongint u=-1

for n as ulongint=u-50 to u-1
    dim as ulongint z=rnd*u
    print n; " mod ";z;"  =  "
    print "fb",,"string"
    print n mod z,mod_(str(n),str(z))
    print
next
sleep

 
Regarding the previuos code, the mod method preserves some granularity.

Code: Select all

 

Dim As Longint lo = 20, hi = 147, R, ans
Dim As Double x
   
R = 7547257594037927936

ans =  R Mod (hi - lo + 1)  + 20 ' dodicat
Print ans
x = R/2^64*(hi - lo + 1) - 0.5 + 20
ans = R/2^64*(hi - lo + 1) - 0.5 + 20  ' dafhi
Print ans, x

print
print
print "Mod","Double"
for z as longint=R-5 to R + 130
 ans =  z Mod (hi - lo + 1)  + lo :print ans ,"   "; 
 ans = z/2^64*(hi - lo + 1) - 0.5 + lo :print ans
 print 
 
 next

Sleep 
The question I ask is not the integrity of mod but:
Both methods put the output in the range uniformly as previously tested.
Which method is more random.
(Tested 64 and 32 bit official fb)
MrSwiss
Posts: 3910
Joined: Jun 02, 2013 9:27
Location: Switzerland

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

Post by MrSwiss »

dodicat wrote:Both methods put the output in the range uniformly as previously tested. Which method is more random.
Not shure to understand properly.
The randomness should come from the gererator (used), not from the ranging
methods process ... ? (as long, as the generated value itself, is used 'as-is')
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] »

@MrSwiss

I wrote a code snippet and tested it with various values. What surprised me was R>=2^32-1. I would have thought R>2^32-1.

Obviously, I did something wrong but I used 'Quick Run' on the snippet and am not able to see what I did wrong.

Now, I go to a lot of effort when posting on this forum to make sure that what I say or code is accurate. The following statement is totally out of order.
This sort of statement is the reason, that I don't trust code of yours.
You remind me of some people who say "You are always doing ....." or "You never do ...." when referring to an event which may never have happened before.
dodicat wrote:mod seems OK for > 2^32-1
Why do you reckon it is not?
That is how you should have responded.

I should tell you that I have just reported your statement to admin. I will say no more until I get a response.
MrSwiss
Posts: 3910
Joined: Jun 02, 2013 9:27
Location: Switzerland

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

Post by MrSwiss »

@deltarho, this must be a joke of yours.

Because, I stayed with it (even if 'sworn in'), by any court of law ...
I simply can't help it, if you can't take the truth.
That's what it is, no more, no less.
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] »

No further comments until admin get back to me.
MrSwiss
Posts: 3910
Joined: Jun 02, 2013 9:27
Location: Switzerland

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

Post by MrSwiss »

deltarho[1859] wrote:No further comments until admin get back to me.
Fine with me. No problem.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

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

Post by dodicat »

Anyway
My question was regarding the randomness.
for simplicity the out values are 0 to 10
Using Dafhi's scenario, it is equivalent to a mapping of all of 0 to 18446744073709551615 into 0 to 10.
So as a number ranges from 0 to 18446744073709551615 another number ranges from 0 to 10.
For many values of a number starting at 0, dafhi's output remains the same until the number passes the first 1/10 th. part of the ulongint range, then the 2/10 th. part e.t.c.

With the mod method every number has a different output as it ranges upwards.
0 to 10 then 0 to 10 ... ....
But because the outputs are 0 to 10 by each method they seem similar, but they couldn't be more different.
Which method produces the best randomness?

Roughly:

Code: Select all



function map(a as double,b as double,x as double,c as double,d as double) as double
    return ((d)-(c))*((x)-(a))/((b)-(a))+(c)
end function


dim as ulongint lo=0,hi=10,y,y2,y3
dim as ulongint R=7547257594037927936
y=(R/2^64)*(hi-lo+1)-.5 +lo
print y;" by dafhi"
y= R Mod (hi - lo + 1)  + lo
print y; " by mod"
y= map(0,2^64-1,R,lo,hi)
print y;" by mapping"
print
print "press a key"
sleep
print

r=0

print "R";tab(30);"dafhi";tab(40);"map";tab(50);"mod"
do
    y=(R/2^64)*(hi-lo+1)-.5 +lo
    y2=map(0,2^64-1,R,lo,hi)
    y3=R Mod (hi - lo + 1)  + lo
    
     print r;tab(30);y;tab(40);y2;tab(50);y3
  r+=  18446744073709551615\10 
    loop until r> 18446744073709551610 or y=10
    sleep

 
  
The mod method happens to produce the same outputs, I wasn't really expecting that.
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

Formula wise the only difference between dafhi and mapping is the divisor being '2^64' and '2^64-1' which will not be seen computationally.

Expanding 'function map' the ordering of computation is different between dafhi and mapping and the difference of 'R' compared with 'R-a' will be negligible for most cases. So, I have dropped mapping and concentrated on dafhi and Mod.

I have changed the divisions from hard-code to '\(hi-lo) and the loop limit from hard-code to 'y = hi'.

Code: Select all

Dim As Ulongint lo=50,hi=100,y,y2,y3
Dim As Ulongint R=7547257594037927936
y=(R/2^64)*(hi-lo+1)-.5 +lo
Print y;" by dafhi"
y= R Mod (hi - lo + 1)  + lo
Print y; " by Mod"
Print
Print "press a key"
Sleep
Print
 
r=0
 
Print "R";Tab(30);"dafhi";Tab(50);"Mod"
Do
  y=(R/2^64)*(hi-lo+1)-.5 +lo
  y3=R Mod (hi - lo + 1)  + lo
  Print r;Tab(30);y;Tab(50);y3
  r+=  18446744073709551615\(hi-lo)            ' \(hi-lo)
Loop Until r> 18446744073709551610 Or y=hi     ' y=hi
Sleep
With 'lo = 0' and 'hi = 10' we get as above.

However, with 'lo = 50' and 'hi = 100' dafhi has a clear mapping correspondence but Mod is all over the place. Small ranges hide the difference. A 'lo = 20' and 'hi = 50' sees Mod locking in at 20 for all r which makes me suspicious of my alterations.

Code: Select all

70 by dafhi
58 by Mod
 
press a key
 
R                            dafhi               Mod
0                            50                  50
368934881474191032           51                  65
737869762948382064           52                  80
1106804644422573096          53                  95
1475739525896764128          54                  59
1844674407370955160          55                  74
2213609288845146192          56                  89
2582544170319337224          57                  53
2951479051793528256          58                  68
3320413933267719288          59                  83
3689348814741910320          60                  98
4058283696216101352          61                  62
4427218577690292384          62                  77
4796153459164483416          63                  92
5165088340638674448          64                  56
5534023222112865480          65                  71
5902958103587056512          66                  86
6271892985061247544          67                  50
6640827866535438576          68                  65
7009762748009629608          69                  80
7378697629483820640          70                  95
7747632510958011672          71                  59
8116567392432202704          72                  74
8485502273906393736          73                  89
8854437155380584768          74                  53
9223372036854775800          75                  68
9592306918328966832          76                  83
9961241799803157864          77                  98
10330176681277348896         78                  62
10699111562751539928         79                  77
11068046444225730960         80                  92
11436981325699921992         81                  56
11805916207174113024         82                  71
12174851088648304056         83                  86
12543785970122495088         84                  50
12912720851596686120         85                  65
13281655733070877152         86                  80
13650590614545068184         87                  95
14019525496019259216         88                  59
14388460377493450248         89                  74
14757395258967641280         90                  89
15126330140441832312         91                  53
15495265021916023344         92                  68
15864199903390214376         93                  83
16233134784864405408         94                  98
16602069666338596440         95                  62
16971004547812787472         96                  77
17339939429286978504         97                  92
17708874310761169536         98                  56
18077809192235360568         99                  71
18446744073709551600         100                 86
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

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

Post by dodicat »

Thanks deltarho.
You would expect mod to be all over the place in almost any range and ulongint split accordingly, I was surprised to see it otherwise in 0 to 10.
Also range 20,30 and 40 and some others see mod a constant regurgitation.
It is all a bit like dafhi's method being akin to a large shell and the mod method a gatling gun.
Both causing equal carnage in different ways.
Anyway, I'll leave it at that for now.
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

It is amazing how one piece of code tells us that we are on a winner and another tells us that all hell has broken loose. I think that we can agree that dafhi's method is more robust. I had already changed the opening post back to dafhi's method. I should tell you that there are issues with linear mapping in general and it is often the case of employing damage control.

I will also leave it a that.
dafhi
Posts: 1640
Joined: Jun 04, 2005 9:51

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

Post by dafhi »

i'm testing the mod method on a ubyte LCG. just need "a few hours" debugging my test.

i kind of already know the result but .. gotta code it

...

Boom.

Code: Select all

function LCG as ubyte:  const mul = 23, add = 3
  static as ushort state:  state = mul * state + add
  return state shr 8
End Function

function range(hi as single=1, lo as single=0)as ubyte
  return lcg mod (hi+1-lo) + lo
End Function


dim as long bins(255)

for i as long = 1 to 256 * 100
  bins( lcg ) += 1
Next

var vlo = bins(0), vhi = vlo
for i as long = 1 to 255
  if vlo > bins(i) then vlo = bins(i)
  if vhi < bins(i) then vhi = bins(i)
  bins(i) = 0
Next: bins(0) = 0

? "low count", "high count";  "  ...  normal range"
? vlo, vhi
?
?
? "halves sum, from range midpoint"
?
? "low", "high", "midpoint"
for range_hi as long = 10 to 30 step 2
 
  for i as long = 1 to 256 * 19
    bins( range( range_hi ) ) += 1
  Next

  dim as ulong m = (range_hi) \ 2
  for i as long = 1 to m
    bins(0) += bins(i)
    bins(range_hi) += bins(range_hi-i)
  Next
 
  ? bins(0), bins(range_hi), m
 
  for i as long = 0 to range_hi
    bins(i)=0
  next
 
next:  sleep
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

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

Post by dodicat »

Hi dafhi.
Yes, it makes sense to use ubyte.
18446744073709551615 is really beyond comprehension for testing.
But I see you have range hi to lo, was that a deliberate trap, it caught me out at first anyhow (but that is not difficult)
Here is the crt rand(), the top half of short on this box.
There seems to be little difference in using the double method and the mod method, only speed.
We can only assume that rand() has been tried and tested for randomness.

Code: Select all

#include "crt.bi"

type range
    as double _max,_min,sd
    as long _maxi,_mini
end type

function mean(a() as double,R as range) as double
    R=type(-1e10,1e10,0,0,0)
    dim as long lb=lbound(a),ub=ubound(a)
    dim as double acc,m
    for n as long=lb to ub
        acc+=a(n)
        if R._max<a(n) then R._max=a(n):R._maxi=n
        if R._min>a(n) then R._min=a(n):R._mini=n
    next
    m=acc/(ub-lb+1)
    acc=0
    for n as long=lb to ub
        acc+=(a(n)-m)*(a(n)-m)
    next
    R.sd =sqr(acc/(ub-lb))
    return m
end function

function chisquared(a() as double,R as range,byref m as double) as double
     m=mean(a(),R)
    dim as double acc
    for n as long=lbound(a) to ubound(a)
        acc+=(a(n)-m)^2/m
    next
    return acc
end function

function rangedbl(lo as long=0, hi as long=1) as long
  return (rand()/rand_max)*(hi+1-lo) -.5 + lo
End Function

function rangemod(lo as long=0, hi as long=1) as long
  return (rand() mod (hi-lo +1))+ lo
End Function

srand(timer)
'================   TEST ============
dim as long lim=5000000
dim as long lo=0
dim as long hi=1

redim as double a(1 to lim),b(1 to lim)
for n as long=1 to lim
    a(n)=(rangemod(lo,hi))
    b(n)=(rangedbl(lo,hi))
next

for n as long=lbound(a) to 20
    print n, a(n),b(n)
next
print "..."


print

dim as range R
dim as double m
var c=chisquared(a(),R,m)
print "Mod"
print "Mean ";m;"  Standard dev. ";R.sd;" Min ";R._min;" Max ";R._max
print "chi squared distance ";c
print
print
print "double"
c=chisquared(b(),R,m)
print "Mean ";m;"  Standard dev. ";R.sd;" Min ";R._min;" Max ";R._max
print "chi squared distance ";c
print
sleep


  
Post Reply