posted February 05, 2007 09:57 AM            

#COMPILE EXE
#DIM ALL
'
MACRO initializeRandomGen       'You choose these 3 parameters at startup as you like, following the rules below:
    soph  = &h68099D27          'Choose from list of Sophie Germain primes (or even find your own!). Maximum start value = &h7fffffff
    carry = &h7eaaaacc          'Maximum start value = &h7ffffffe and always choose it to be less than soph (and not less than 3)
    mwc   = &he6690adb          'You choose any dword 3 to &hffffffff (4294967295)
END MACRO
' here's a little list of a few more soph primes:
' &h7C6BAD4C;&h6847AFAF;&h6844E021;&h686CDDA7;&h68661E4F;&h7C6BCEFD;&h7C6AAD74
' Each one produces its own unique and independent sequence of about 2^62 random dwords, and there are thousands of sophs out there.
' If x * 2^32 - 1 is prime and x * 2^31 - 1 is prime, then x is good to go as a soph for this algorithm.
'
MACRO getRandomDword         'here is the algorithm engine for the random numbers
    tmpProd = soph * mwc     'soph size importance is apparent here. Above &h7fffffff it overloads the temp product.
    tmpProd = tmpProd + carry
    mwc     = @tpPtr         'here is your dword random number. 1 dword (4 bytes) per call
    carry   = @tpPtr[1]      'carry gets saved too for the next random call
END MACRO                    'you can also do this as a MACRO FUNCTION getRandomDword: END MACRO = mwc
'
GLOBAL soph, carry AS LONG      'make these global so you can use the random generator thru all your program.
GLOBAL tmpProd AS QUAD
GLOBAL tpPtr AS LONG PTR
GLOBAL mwc AS DWORD
'                               'PB win 8
FUNCTION PBMAIN () AS LONG      'this is a Multiply with Carry algorithm random # generator. It has passed all randomness tests so far.
                                'Below we create a binary file of random dwords. Period is approx 2^62 dwords
    LOCAL ii, i2 AS LONG
    LOCAL trsPtr AS LONG PTR
    LOCAL testRandString AS STRING
'
    OPEN "rndMwcStandardPBcode.dat" FOR BINARY AS #1
'
    testRandString = SPACE$(400000)
    tpPtr = VARPTR(tmpProd)
'
    initializeRandomGen     'set initial parameters
'
    FOR i2 = 1 TO 100       '100 loops creates a 40MB file
       trsPtr = STRPTR(testRandString)
       FOR ii = 1 TO 100000
          getRandomDword    'mwc (stands for Multiply with Carry) is produced here which is your next random number
          @trsPtr = mwc     'write it to a test-string
          INCR trsPtr       'test-string pointer
       NEXT
       PUT #1, , testRandString
    NEXT
     ? "k"
END FUNCTION

posted February 05, 2007 11:10 AM            

I am wary of such data, partly but not exclusively because
I have found no evidence that estimates of (for example) pi, based
on Random.org data converge any more closely or quickly to theoretically
"expected" values than do estimates based on pseudo-RNG's

posted February 05, 2007 08:37 PM            

#COMPILE EXE
#DIM ALL
                           'PB ver. 8
FUNCTION PBMAIN () AS LONG 'This program first creates a 4GB file across the entire PB rnd generator period, 1 byte per RND,
   'using 5777766666 as the seed. Then it will re-seed and try to match the first 8, then 65536 more bytes in a row, and if all
   'match, it verifies that different seeds produce the same 4GB series, but shifted to another point in the series. Run
   'enough different seeds, and it becomes assured that all seeds produce the same 4GB series, with a shift.
'
    LOCAL ii, i2, i3 AS LONG
    LOCAL eRnd AS EXT, t1 AS DOUBLE
    LOCAL erPtr, trsPtr, srsPtr AS BYTE PTR
    LOCAL testRandString AS STRING * 65536
    LOCAL shortRandString AS STRING * 8
    LOCAL lineo AS STRING
'
    RANDOMIZE 5777766666           'any valid number here can be used
    OPEN "d:\rndBytePeriod.dat" FOR BINARY AS #1
'
    erPtr  = VARPTR(eRnd)          'ptr to the extended precision rand numb
'    t1 = TIMER                    'REMOVED. just test code 
'   GOTO reseed                    'once generated, you can go directly to the re-seed code to test further random seeds,
                                   'as it takes several minutes to initially create the full period 4GB file.
    FOR ii = 1 TO 65536
    trsPtr = VARPTR(testRandString)
       FOR i2 = 1 TO 65536
'
          eRnd    = RND
          @trsPtr = @erPtr[6]        'get just let's say the 7th byte of the EXT random number since it's quite a random 1,
                                     'and write it to testRandString
          INCR trsPtr
       NEXT
       PUT #1, , testRandString
    NEXT
'
   '-----------------------------------------------------------------------------------------------------------
    reseed:                          'here we now test to see if we find the same sequences from any seed
'
    RANDOMIZE
    SEEK #1, 1                       'ADDED. Bug hopefully fixed
    srsPtr = VARPTR(shortRandString) 'create random 8 bytes to search for
'
       FOR i2 = 1 TO 8               'you can search for longer than 8 byte strings, but 10^18 or so odds with
                                     'repeatability is beyond reasonable doubt, imho.
          eRnd    = RND
          @srsPtr = @erPtr[6]        'get just the 7th byte of the EXT random number again to check for repeats
'
          INCR srsPtr
       NEXT
       DO
       GET$ #1, &h1000000, lineo
       IF INSTR(lineo, shortRandString) THEN
           ? "Found it! The odds of that happening by chance alone are over" & $CRLF & _
             "1 in a billion billion--about the same as my state lottery     " & $CRLF & _
             "Found at position " & FORMAT$(SEEK(#1) - &h1000000 + INSTR(lineo, shortRandString))
           ? "I will now generate 65536 more randoms in sequence and see if they all match.."
'
           trsPtr = VARPTR(testRandString)
           SEEK #1, SEEK(#1) - &h1000000 + INSTR(lineo, shortRandString) + 7
           GET #1, , testRandString
              FOR i3 = 1 TO 65536
'
                 eRnd    = RND
                 IF @trsPtr <> @erPtr[6] THEN ? "Not the same! But, may have reached end of file"
'
                 INCR trsPtr
              NEXT
              ? "All matched! Both sequences are the same, no question now."
              EXIT FUNCTION
       END IF
       LOOP UNTIL EOF(1)
'
       ? "Didn't find the sequence! This is probably because it is on a boundary between string reads, which for speed" & _
       "reasons and rarity of occurrence, I decided not to account for."
'
END FUNCTION

posted February 06, 2007 01:43 AM            

What did you intend doing with t1?

The app finished without saying a word to me.

Anyway, just immediate reactions - I'll read it better later.

Added: Ran it again and it fell short again but this time I got "Didn't find the sequence! This ..."

posted February 06, 2007 05:03 AM            

Using seed values from 3 to 1999, I compared PB'sRND(), John Gleason's and Blum Blum Shub.
As before we are comparing the 'coin flip' scores with a binomial distribution and applying
a Chi-squared on each measurable category. If the Chi exceeds 35 (on 30 degrees of freedom),
then the seed is rejected.

Each seed is tested with 26*1000*30 = 780,000 flips

PB RND() 26 seeds
John Gleason 28 seeds
Blum Blum Shub 38 seeds

I also looked at a 10Mb random bit file from
http://www.random.org

as Emil suggested, but to carry out a comparable test would involve
780,000 * 2000/8 about 200 megs of downloaded data. So I just looked at 1 file
and found the profiles looked very similar in 'noisiness' to pseudo-random
numbers.

My code looks very messy after all that experimenting but if any one would like
to see it. I can post it here, in Powerbasic this time! Script is 60 times slower.

posted February 06, 2007 07:52 AM            

The timer variable is residual test code I forgot to remove. oof.

The second problem, however, is big trouble because it refutes my
proof. When I ran the code, it always found the sequences every
time with dozens of different seeds. Let me ask this: when you run the
second section of the program after generating the 4GB file, have you
used the GOTO to check other seeds? Is the program failing to find the
sequence immediately or is it clearly searching thru the 4GB file and
still failing to find it?

ADDED: Geez sorry all, I found the bug. I'll add a line to reset the
4GB file to its beginning before checking new seeds.

posted February 06, 2007 08:53 AM            

It's a HUGE challenge to get a good physical rng without applying
some randomizer like AES to the output. But compression of the output
(without password of course) I do consider "fair means."

posted February 06, 2007 10:56 AM            

One improvement that could be made is to lump the se loop and
the g loop into one, then reduce the Chi degrees of freedom to
1 instead of 30.

I have configured it to use your algorithm John, but the others
are there, commented out.

#COMPILE EXE
#DIM ALL
'
' Charles E V Pegge
' 6 February 2007
'
'
' Testing Random Number Generators
'
' This will take several minutes to run and produces
' a smallish file called 'results.txt' with the selected
' seeds and their profiles
'
' use as you please
'
'
'
'
'
'
' John Gleason's Random Number Generator Macros
' http://www.powerbasic.com/support/forums/Forum12/HTML/004081-8.html 
' 5 February 2007
'
MACRO initializeRandomGen 'You choose these 3 parameters at startup as you like, following the rules below:
 soph  = &h68099D27 'Choose from list of Sophie Germain primes (or even find your own!). Maximum start value = &h7fffffff
 carry = &h7eaaaacc 'Maximum start value = &h7ffffffe and always choose it to be less than soph (and not less than 3)
 mwc   = seed ' &he6690adb 'You choose any dword 3 to &hffffffff (4294967295)
END MACRO
'
' here's a little list of a few more soph primes:
' &h7C6BAD4C;&h6847AFAF;&h6844E021;&h686CDDA7;&h68661E4F;&h7C6BCEFD;&h7C6AAD74
' Each one produces its own unique and independent sequence of about 2^62 random dwords, and there are thousands of sophs out there.
' If x * 2^32 - 1 is prime and x * 2^31 - 1 is prime, then x is good to go as a soph for this algorithm.
'
MACRO getRandomDword 'here is the algorithm engine for the random numbers
 tmpProd = soph * mwc 'soph size importance is apparent here. Above &h7fffffff it overloads the temp product.
 tmpProd = tmpProd + carry
 mwc     = @tpPtr 'here is your dword random number. 1 dword (4 bytes) per call
 carry   = @tpPtr[1] 'carry gets saved too for the next random call
END MACRO 'you can also do this as a MACRO FUNCTION getRandomDword: END MACRO = mwc
'
GLOBAL soph, carry AS LONG 'make these global so you can use the random generator thru all your program.
GLOBAL tmpProd AS QUAD
GLOBAL tpPtr AS LONG PTR
GLOBAL mwc AS DWORD
'
'------
'
GLOBAL ranbuf AS STRING
GLOBAL ranbufp AS BYTE PTR
GLOBAL lenranbuf AS LONG
GLOBAL byto AS LONG
GLOBAL byti AS LONG
GLOBAL biti AS LONG
'
' Accessing a string of random bits, loaded from a file
FUNCTION ranbit() AS LONG
 IF biti=0 THEN byti=@ranbufp[byto]: byto=(byto+1) MOD lenranbuf
 FUNCTION=BIT(byti,biti)
 biti=(biti+1) MOD 8
END FUNCTION
'
' used for handling large factorials
FUNCTION lfactorial(BYVAL b AS LONG) AS DOUBLE
IF b < 2 THEN EXIT FUNCTION
DIM a AS DOUBLE
a=0
DO
 a=a+LOG(b): DECR b: IF b <= 1 THEN EXIT LOOP
LOOP
FUNCTION=a
END FUNCTION
'
'
FUNCTION PBMAIN () AS LONG
'
'
' For John's algorithm
tpPtr = VARPTR(tmpProd)
'
'// chi squared check main routine
'
'
 DIM x AS DOUBLE,g AS DOUBLE, sdf AS DOUBLE, se AS DOUBLE, p AS DOUBLE
 DIM seed AS DOUBLE, seede AS DOUBLE, seedcount AS DOUBLE
 DIM j AS DOUBLE, i AS DOUBLE, t AS DOUBLE, m AS DOUBLE
 DIM d(100) AS DOUBLE, ch(100) AS DOUBLE, bn(100) AS DOUBLE
 DIM pp AS DOUBLE, c AS DOUBLE, sc AS DOUBLE, mx AS DOUBLE
 DIM n AS DOUBLE, r AS DOUBLE, q AS DOUBLE,xx AS DOUBLE, xi AS DOUBLE
 DIM rk AS QUAD, rm AS QUAD, ra AS QUAD
 DIM rmb AS QUAD,rab AS QUAD
'
'
 OPEN "results.txt" FOR OUTPUT AS 1
'
 ' for sample from http://www.random.org 
 ' OPEN "ran1.bin" FOR BINARY AS 2
 ' GET$#2,LOF(2),ranbuf
 ' CLOSE(2)
 ' ranbufp=STRPTR(ranbuf)
 ' lenranbuf=LEN(ranbuf)
 ' PRINT#1,"Random bit file length "+STR$(lenranbuf)+" bytes"
 ' PRINT#1,""
'
 ' our run parameters
 x=26:g=1e3: sdf=30: se=sdf: p=.5: seed=3: seede=2000
'
 j=g: m=x+1: pp=0: mx=0
 n=x: r=0: q=1-p: xx=0: xi=0
'
 FOR i=0 TO m:d(i)=0:ch(i)=0:bn(i)=0:NEXT
'
'
PRINT #1," Chi Squared Variance Probability Distribution for p=" +STR$(p)+" with "+STR$(g)+" iterate loops of "+STR$(x)+" in "+STR$(sdf)+" sessions"
PRINT #1," Comparing with binomial distribution for ("+STR$(p)+"+"+STR$(1-p)+")^"+STR$(x)
PRINT #1," "+STR$(se)+" degrees of freedom for each of "+STR$(x)+" categories.
PRINT #1," Randomizer seed="+STR$(seed)
'
' prepare binomial table to compare with
r=0
DO
  bn(r)=EXP(lfactorial(n)-lfactorial(r)-lfactorial(n-r) ) *p^(n-r)*q^r*g
  INCR r: IF r>n THEN EXIT LOOP
LOOP
'
'
' test a seed
DO
'
 'Sowing the seed
'
 ' PB internal
 RANDOMIZE(seed)
'
 ' John Gleason
 ' macro modified to read 'seed' into mwc
 initializeRandomGen
'
 ' BLUM BLUM SHUB
 ' 2 primes both with mod 4 equal to 3
 rmb=32027*64007: rab=seed
'
 ' this one is too big for a QUAD
 rm=2^48-1: rk=31167285: ra=seed
'
'
 DO
  ' run a batch of se sessions
  j=g
  ' run 1000 sets of throws
  DO
   i=x: t=0
   DO
   ' 1 set of x throws
'
    ' Various methods of getting a random 0/1 to add to t
'
    ' ra=(ra*rk) mod rm: if ra and 2 then incr t
'
    ' Blum Blum Shub RNG
    '  rab=(rab*rab) MOD rmb
    ' to use the bits in the lowest order byte
    '  IF biti=0 THEN getRandomDword
    '  IF BIT(rab,biti) THEN INCR t
    '  biti=(biti+1) MOD 8
'
    ' normal PB method
    ' IF RND(1) >= p THEN INCR t
'
    ' John Gleason RNG
    ' to use the bits in the lowest order byte
     IF biti=0 THEN getRandomDword
     IF BIT(mwc,biti) THEN INCR t
     biti=(biti+1) MOD 8
'
    'from random number file
    'if ranbit() then incr t
'
    DECR i:IF i <= 0 THEN EXIT LOOP
   LOOP
   ' use the score to index the right 'bin' then increment it
   INCR d(t+1): DECR j: IF j <= 0 THEN EXIT LOOP
  LOOP
'
  ' after a run of g sets of x throws, the chi data is accummulated
  i=m: c=1: r=0
  DO
   pp=bn(r) ' expected result according to binomial distribution
   xi=( (d(c)-pp)^2/pp ) + ( ( (g-d(c))-(g-pp) )^2 / (g-pp) )  '// chi terms
   IF pp >= 2 THEN xx=xx+xi: ch(m-i)=ch(m-i)+xi ' exclude very low Chi
'
   ' reject this seed if a chi value exceeds 35
   ' this applies to se=30 for 30 degrees of freedom
   ' please refer to a Chi-squared chart.
   IF ch(m-i) > 35 THEN GOTO another_seeding
'
   d(c)=0: INCR r: INCR c: DECR i: IF i <= 0 THEN EXIT LOOP
  LOOP
  DECR se: IF se <= 0  THEN EXIT LOOP
 LOOP
'
 ' work out a scale for the bar chart
 i=m: mx=0
 DO
  mx=MAX(mx,ch(i)): DECR i: IF i < 0 THEN EXIT LOOP
 LOOP
'
 ' print the results
 i=m: sc=60/mx
 PRINT#1,""
 PRINT #1," New seed "+STR$(seed): INCR seedcount
 PRINT#1,""
 DO
  PRINT #1," |"+REPEAT$(sc*ch(m-i),"_")+"  #"+USING$("##",m-i)+" ChiSq ="+USING$("##.#",ch(m-i))
  DECR i: IF i <= 0 THEN EXIT LOOP
 LOOP
 PRINT #1,""
 PRINT #1," total chiSq="+STR$(xx)+"  at "+STR$(x*sdf)+" degrees of freedom"
'
'
 ' prepare for the next seed to process
 another_seeding:
 IF seed > seede THEN EXIT LOOP
'
 'clear the data arrays for another seed test
 i=m
 DO
  ch(i)=0:d(i)=0: DECR i:IF i < 0 THEN EXIT LOOP
 LOOP
 INCR seed
 se=sdf
 xx=0
'
' back to the top
LOOP
'
PRINT#1,""
PRINT#1,""
PRINT#1,"Total seed count: "+STR$(seedcount)
'
CLOSE(1)
'
END FUNCTION

------------------
www.pegge.net

posted February 07, 2007 01:27 AM            

quote:

---

I should be looking at the actual distribution compared with the expected one much along the lines that Charles has been doing with the binomial and Chi-squared.

---

 1% 1.0036%
 2% 2.0028%
 4% 4.0023%
 8% 7.9936%
16% 15.9935%
32% 32.0031%
50% 49.9938%
50% 50.0062%
32% 32.0099%
16% 15.9789%
 8% 7.9988%
 4% 4.0063%
 2% 2.0057%
 1% 0.9993%

Charles, I haven't had a chance to run your code yet. What conclusions are you drawing so far?

posted February 07, 2007 06:39 AM            

1 Most seeds generate pseudo random numbers which show patterns of excessive deviation
from the expected norm.

2 Seeds can be selected that avoid such excesses, within a specified amount of random numbers.

3 This could be useful if you need your rands to be bland.

4 Random is a theoretical concept, rarely to be found in nature, since everything is
caused by something else. So selecting a theoretically random set numbers is not
a sin.

posted February 07, 2007 08:27 AM            

quote:

---

' reject this seed if a chi value exceeds 35
' this applies to se=30 for 30 degrees of freedom
' please refer to a Chi-squared chart.

---

posted February 07, 2007 09:25 AM            

John, your program works fine (after I added a couple of WAITKEY$
so that I could read the printed output) and I am impressed. I have
one question: Do the data in your 4GB file contain all the EXT numbers
produced by the RND function -- so that, if I extract 8 bytes at a time
& convert them to EXT's, I can study them further?

I would be more expansive, but am pressed for time right now.

posted February 07, 2007 10:48 AM            

I have just looked at at proper Chi table and I agree that my
Chi limit was set too low. If you change this figure from 35 to 43.8,
this will represent a .05 confidence limit. You'll certainly
get more seeds! source:
http://www.statsoft.com/textbook/sttable.html

I would like to understand the basis of Chi better, going back
to first principles would be nice. But for now I am assuming it is a
form of variance with a standard deviation of 1 (@ 1df).

posted February 07, 2007 01:17 PM            

quote:

---

1 Most seeds generate pseudo random numbers which show patterns of excessive deviation
from the expected norm.

---

xi=( (d(c)-pp)^2/pp ) + ( ( (g-d(c))-(g-pp) )^2 / (g-pp) )  '// chi terms

So, xi = ( (d(c)-pp)^2/pp ) + ( (d(c)-pp)^2/ (g-pp) )
       = (d(c)-pp)^2 * ( 1/pp + 1/(g-pp) )

posted February 07, 2007 02:10 PM            

With a little modification, you can get each entire 10-byte EXT number written
to the file so you can look at it. Here's a little stand-alone modification
that can make a really big file of the full period of EXT rand numbers if wanted,
but it's commented out to limit it initially to 4GB:

#COMPILE EXE
#DIM ALL
'
FUNCTION PBMAIN () AS LONG 'This can create a huge 40GB file of EXT randoms
'
    LOCAL ii, i2 AS LONG
    LOCAL eRnd AS EXT
    LOCAL erPtr, trsPtr         AS EXT PTR   '<<CHANGED to extended ptr
    LOCAL testRandString AS STRING * 655360  '<<CHANGED ten times bigger than before
    RANDOMIZE 5777766666           'any valid number here can be used, or blank/TIMER
    OPEN "c:\rndEXTperiod.dat" FOR BINARY AS #1 '<<CHANGED
'
    erPtr  = VARPTR(eRnd)          'ptr to the extended precision rand numb
'
    FOR ii = 1 TO 6553'6           'make this the FULL 40GB PERIOD by uncommenting the '6
    trsPtr = VARPTR(testRandString)
       FOR i2 = 1 TO 65536
'
          eRnd    = RND
          @trsPtr = @erPtr        '<<CHANGED write the whole 10 yards to testRandString
          INCR trsPtr
       NEXT
       PUT #1, , testRandString
    NEXT
'
END FUNCTION