RND breadth

General FreeBASIC programming questions.
albert
Posts: 5953
Joined: Sep 28, 2006 2:41
Location: California, USA

Re: RND breadth

Postby albert » Feb 03, 2019 18:35

The problem is in the RND Function... it only generates upper values... > .0010
Here's a test to see: it only generates out of 1,000,000 loops , it only generates 1,000 below .0010

Code: Select all


screen 19

dim as double l0 = .00000 , h0 = .0001
dim as double l1 = .00011 , h1 = .0002
dim as double l2 = .00021 , h2 = .0003
dim as double l3 = .00031 , h3 = .0004
dim as double l4 = .00041 , h4 = .0005
dim as double l5 = .00051 , h5 = .0006
dim as double l6 = .00061 , h6 = .0007
dim as double l7 = .00071 , h7 = .0008
dim as double l8 = .00081 , h8 = .0009
dim as double l9 = .00091 , h9 = .0010

dim as integer array(0 to 9)

dim as double n = 0
for a as longint = 1 to 1000000
   
    n = rnd
   
    select case n
    case l0 to h0 : array(0)+=1
    case l1 to h1 : array(1)+=1
    case l2 to h2 : array(2)+=1
    case l3 to h3 : array(3)+=1
    case l4 to h4 : array(4)+=1
    case l5 to h5 : array(5)+=1
    case l6 to h6 : array(6)+=1
    case l7 to h7 : array(7)+=1
    case l8 to h8 : array(8)+=1
    case l9 to h9 : array(9)+=1
    end select
   
next

dim as longint total = 0
for a as longint = lbound(array) to ubound(array) step 1
    print array(a)
    total+=array(a)
next
print "total = " ; total

sleep
end

dodicat
Posts: 6809
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: RND breadth

Postby dodicat » Feb 03, 2019 19:18

.001 is one thousandth of the way from 0 to 1.
So each leg you expect 100, and you have 10 legs up to .001.(as you see printed)
You expect about 1000 below .001 because you have 1000 of these intervals between 0 and one, 1000*1000 = 1,000,000
Here is the independent random generator and I have closed up your gaps so the intervals butt on to each other.

Code: Select all

 
 
'screen 19
' ulongint generator

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 shared as rand z
init(z)
'========
randomize ,0

screen 19

dim as double l0 = .00000 , h0 = .0001
dim as double l1 = .0001 , h1 = .0002
dim as double l2 = .0002 , h2 = .0003
dim as double l3 = .0003 , h3 = .0004
dim as double l4 = .0004 , h4 = .0005
dim as double l5 = .0005 , h5 = .0006
dim as double l6 = .0006 , h6 = .0007
dim as double l7 = .0007 , h7 = .0008
dim as double l8 = .0008 , h8 = .0009
dim as double l9 = .0009 , h9 = .0010

dim as integer array(0 to 9)

dim as double n = 0
for a as longint = 1 to 1000000
   
    n = randouble(z)
   
    select case n
    case l0 to h0 : array(0)+=1
    case l1 to h1 : array(1)+=1
    case l2 to h2 : array(2)+=1
    case l3 to h3 : array(3)+=1
    case l4 to h4 : array(4)+=1
    case l5 to h5 : array(5)+=1
    case l6 to h6 : array(6)+=1
    case l7 to h7 : array(7)+=1
    case l8 to h8 : array(8)+=1
    case l9 to h9 : array(9)+=1
    end select
   
next

dim as longint total = 0
for a as longint = lbound(array) to ubound(array) step 1
    print array(a)
    total+=array(a)
next
print "total = " ; total

sleep
end
 
albert
Posts: 5953
Joined: Sep 28, 2006 2:41
Location: California, USA

Re: RND breadth

Postby albert » Feb 03, 2019 19:20

The RND stuff is for my program Vari_Cph...

I needed to generate a 64 bit random binary string...
But it only generated upper values ,mostly "1's"

Right( string( 64 , "0" ) + bin( int( rnd * (2^64) ) ) , 64 ) , was generating mostly 1's..

I got it returning a full plethora of bits...With the following code..

Code: Select all


screen 19

do
   
    dim as string*64 str1 = right(string(64,"1"),64) ' + bin(int(rnd*2^64)),64)
   
    for a as longint = 1 to 48
        mid(str1,int(rnd*len(str1))+1,1) = "0"
    next
   
    print str1 , valulng("&B"+str1)

    sleep

loop until inkey = chr(27)

end
   
sancho3
Posts: 358
Joined: Sep 30, 2017 3:22

Re: RND breadth

Postby sancho3 » Feb 03, 2019 19:23

Like Dodicat said, you have isolated a small range of numbers. If you isolate a same size range anywhere else you get similar amounts of hits.
Here I add the next range .00100 to .00200 to your code and we see similar amounts of hits:

Code: Select all

screen 19

dim as double l0 = .00000 , h0 = .0001
dim as double l1 = .00011 , h1 = .0002
dim as double l2 = .00021 , h2 = .0003
dim as double l3 = .00031 , h3 = .0004
dim as double l4 = .00041 , h4 = .0005
dim as double l5 = .00051 , h5 = .0006
dim as double l6 = .00061 , h6 = .0007
dim as double l7 = .00071 , h7 = .0008
dim as double l8 = .00081 , h8 = .0009
dim as double l9 = .00091 , h9 = .0010

dim As Double ll = .00100, hh = .0020

dim as integer array(0 to 10)

dim as double n = 0
for a as longint = 1 to 1000000
   
    n = rnd
   
    select case n
    case l0 to h0 : array(0)+=1
    case l1 to h1 : array(1)+=1
    case l2 to h2 : array(2)+=1
    case l3 to h3 : array(3)+=1
    case l4 to h4 : array(4)+=1
    case l5 to h5 : array(5)+=1
    case l6 to h6 : array(6)+=1
    case l7 to h7 : array(7)+=1
    case l8 to h8 : array(8)+=1
    case l9 to h9 : array(9)+=1
    Case ll To hh : array(10)+= 1
    end select
   
next

dim as longint total = 0
for a as longint = lbound(array) to ubound(array) step 1
    print array(a)
    total+=array(a)
next
print "total = " ; total
Print "next range: "; array(10)
sleep
end

deltarho[1859]
Posts: 2851
Joined: Jan 02, 2017 0:34
Location: UK

Re: RND breadth

Postby deltarho[1859] » Feb 03, 2019 19:36

Oh dear!

dodicat and sancho3 beat me to it.

Lets look at the first two lines of the Select construct.

Code: Select all

Case l0 To h0 => 0 To .0001
Case l1 To h1 => .00011 To .0002

What happens to the values > .0001 and < .00011 ie a span of .00001. Simple - they don't get counted.

The Mersenne Twister passes a PractRand test to 256GB. If your 'problem', albert, really did exist PractRand would fail long before 256GB.

I am definitely out of this thread now. <smile>
jj2007
Posts: 1961
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: RND breadth

Postby jj2007 » Feb 03, 2019 19:41

albert wrote:But it only generated upper values ,mostly "1's"
If you use a 2^64 range, i.e. 0 ... 18446744073709551615, then the highest number in the first 0.1% of that range is 18446744073709551.615 (but there are probably a dozen different ways to explain the phenomenon - unfortunately it is not very intuitive).
dodicat
Posts: 6809
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: RND breadth

Postby dodicat » Feb 03, 2019 19:57

Albert.
I would create the zeros and ones this way:

Code: Select all

randomize
screen 19

#define range(f,l) Int(Rnd*(((l)+1)-(f)))+(f)

do
    dim as integer countzero,countone

    dim as string str1 = string(64,0)
   
    for a as longint = 0 to len(str1)-1
        var z=range(48,49)
        if z=48 then countzero+=1
        if z=49 then countone+=1
       
        str1[a]=z
    next
   
    print str1 , valulng("&B"+str1)


print "zeros   ";countzero
print "ones    ";countone
 if countone+countzero<>64 then print "Error!"
    sleep

loop until inkey = chr(27)

end
     
counting_pine
Site Admin
Posts: 6245
Joined: Jul 05, 2005 17:32
Location: Manchester, Lancs

Re: RND breadth

Postby counting_pine » Feb 03, 2019 20:31

albert wrote:The RND stuff is for my program Vari_Cph...

I needed to generate a 64 bit random binary string...
But it only generated upper values ,mostly "1's"

I guess the problem you're having is that Rnd can't generate 64 bits of randomness in each number.
It should be able to manage up to 53 bits precision because it's a Double, but in practice, it seems to only return up to 32 bits of randomness.
Either way, a full 64 bits isn't possible.

The easiest way is probably just to generate two 32-bit numbers, and combine them:

Code: Select all

function rnd32() as ulong
    return int(rnd * 2^32)
end function

function rnd64() as ulongint
    return culngint(rnd32()) shl 32 + rnd32()
end function

randomize
print bin(rnd64, 64)
deltarho[1859]
Posts: 2851
Joined: Jan 02, 2017 0:34
Location: UK

Re: RND breadth

Postby deltarho[1859] » Feb 04, 2019 14:33

I have popped back to see what counting_pine posted. <smile>
it seems to only return up to 32 bits of randomness.

That is because the underlying FB generators are only 32 bit.

Since we are only requesting a couple of random numbers speed is not an issue so I would recommend 'Randomize , 5' before calling rnd32(). If our application uses random numbers for other purposes then we can always invoke another Randomize using one of the other PRNGs after calling rnd64().

I actually use something very similar to counting_pine's code.

Code: Select all

#define Get64Bit (Cast( Ulongint, Rnd*(2^32) ) Shl 32) Or Cast( Ulongint, Rnd*(2^32) )

Code: Select all

Randomize , 5
Dim As Ulongint x = Get64Bit
Print Bin(x,64)
Sleep
speedfixer
Posts: 468
Joined: Nov 28, 2012 1:27
Location: California

Re: RND breadth

Postby speedfixer » Jul 20, 2020 8:05

Just stumbled across this.
I know this is old and stale, but I hate to think that Albert never got comfortable with anyone's reason why he did not see what he expected.

Albert:

2^32 is not half of 2^64.

2^63 is half

2^62 is 1/4
2^61 is 1/8
...
...
...
2^32 is VERY small next to 2^64.

If we allow a perfect rnd generator, anything less than 2^32 simply will not be chosen frequently in the 2^64 range.

I hope this helps someone.

david
deltarho[1859]
Posts: 2851
Joined: Jan 02, 2017 0:34
Location: UK

Re: RND breadth

Postby deltarho[1859] » Jul 20, 2020 15:28

2^k * 2^n = 2^(k+n)

So, 2^64 * 1/2 = 2^64 * 2^(-1) = 2^63
MrSwiss
Posts: 3726
Joined: Jun 02, 2013 9:27
Location: Switzerland

Re: RND breadth

Postby MrSwiss » Jul 20, 2020 21:08

Instead of a lot of words, some code to visualize things:

Code: Select all

' ULI_vs_UL_visualized.bas -- (c) 200-07-20, MrSwiss

Width 120, 25
Dim As ULongInt ULI = &hFFFFFFFFFFFFFFFF
Dim As ulong    UL  = &hFFFFFFFF

Print "64 bits vs. 32 bits visualized"
Print "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
Print

Print "64 bits max:  "; Tab(15); Bin(ULI, 64), ULI
ULI = BitReset(ULI, 63)
Print "halve value:  "; Tab(15); Bin(ULI, 64), ULI
ULI = BitReset(ULI, 62)
Print " 4th value:   "; Tab(15); Bin(ULI, 64), ULI
ULI = BitReset(ULI, 61)
Print " 8th value:   "; Tab(15); Bin(ULI, 64), ULI
ULI = BitReset(ULI, 60)
Print "16th value:   "; Tab(15); Bin(ULI, 64), ULI
ULI = BitReset(ULI, 59)
Print "32th value:   "; Tab(15); Bin(ULI, 64), ULI
ULI = BitReset(ULI, 58)
Print "64th value:   "; Tab(15); Bin(ULI, 64), ULI
Print
Print "32 bits max:  "; Tab(47); Bin(UL, 32), UL

Sleep
Result:

Code: Select all

64 bits vs. 32 bits visualized
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

64 bits max:  1111111111111111111111111111111111111111111111111111111111111111      18446744073709551615
halve value:  0111111111111111111111111111111111111111111111111111111111111111      9223372036854775807
 4th value:   0011111111111111111111111111111111111111111111111111111111111111      4611686018427387903
 8th value:   0001111111111111111111111111111111111111111111111111111111111111      2305843009213693951
16th value:   0000111111111111111111111111111111111111111111111111111111111111      1152921504606846975
32th value:   0000011111111111111111111111111111111111111111111111111111111111      576460752303423487
64th value:   0000001111111111111111111111111111111111111111111111111111111111      288230376151711743

32 bits max:                                  11111111111111111111111111111111      4294967295
deltarho[1859]
Posts: 2851
Joined: Jan 02, 2017 0:34
Location: UK

Re: RND breadth

Postby deltarho[1859] » Jul 20, 2020 21:35

For the life of me I cannot see how that will help anyone.

2^k * 2^n = 2^(k+n), says it all.

To keep in the spirit of the last post, whatever that is, just replace all the BitReset() with 'ULI Shr 1'.

Return to “General”

Who is online

Users browsing this forum: Google [Bot] and 8 guests