Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.forth > #26933 > unrolled thread

Multiply With Carry random generators

Back to article view | Back to comp.lang.forth

Contents

  Multiply With Carry random generators albert@spenarnc.xs4all.nl (Albert van der Horst) - 2013-11-25 01:44 +0000
    Re: Multiply With Carry random generators mhx@iae.nl - 2013-11-25 10:40 -0800
    Re: Multiply With Carry random generators C G Montgomery <cgm@physics.utoledo.edu> - 2013-11-25 13:46 -0500

#26933 — Multiply With Carry random generators

Looking into random number generators lately, mhx mentioned the
diehard program.
This was written by Marsaglia and with his linux version of the program he
distributes a file mwcp.ps explaining the Multiply With Carry random generators.

They are simple enough to be typed in by heart whenever you need one,
yet they have respectable periods of 2^64 for 32-bits and 2^128 for
64 bits. Where Marsaglia promotes them they must be pretty good.
The nice thing are the constants. They are small enough to remember.
The multipliers are from wikipedia.
Like the congruentual generators they are simpler in Forth than in other
languages, especially given the double precision wordsets.

\ -------------- untested code ---------------------

WANT [IF]

1 CELLS 4 = [IF]
-178  CONSTANT multiplier  \ 4,294,967,118 32 bits
[ELSE]   \ Ignoring 16 bits.
-742  CONSTANT multiplier   \ 18,446,744,073,709,550,874
[THEN]

CREATE SEED   1 , 1 ,

: RAND SEED 2@   >R multiplier UM*  R> 0 D+ SWAP OVER SEED 2! ;

: RANDOMIZE SEED 2! ;

\ ------------------------------------------------------

I looked at it, and it doesn't seem too shabby.

Groetjes Albert
-- 
Albert van der Horst, UTRECHT,THE NETHERLANDS
Economic growth -- being exponential -- ultimately falters.
albert@spe&ar&c.xs4all.nl &=n http://home.hccnet.nl/a.w.m.van.der.horst

[toc] | [next] | [standalone]

#26946

On Monday, November 25, 2013 2:44:33 AM UTC+1, Albert van der Horst wrote:
 
> Like the congruentual generators they are simpler in Forth than in other
> languages, especially given the double precision wordsets.
> 
> \ -------------- untested code ---------------------
> 
> WANT [IF]
> 
> 1 CELLS 4 = [IF]
> -178  CONSTANT multiplier  \ 4,294,967,118 32 bits
> [ELSE]   \ Ignoring 16 bits.
> -742  CONSTANT multiplier   \ 18,446,744,073,709,550,874
> [THEN] 
> 
> CREATE SEED   1 , 1 ,
> 
> : RAND SEED 2@   >R multiplier UM*  R> 0 D+ SWAP OVER SEED 2! ;
> : RANDOMIZE SEED 2! ;
> \ ------------------------------------------------------
> 
> I looked at it, and it doesn't seem too shabby.

It's one of the best so far. OPERM5 sniffs it out.
It could be worth porting diehard to 64bits for :-)

-marcel
-- -----------------------------
FORTH> diehard
DIEHARD is checking rand-AvdH..

        .---------------------------------------------------------------.
        |           This is the "tough" BIRTHDAY SPACINGS TEST          |
        | Choose 4096 birthdays in a "year" of 2^32 days. Thus each     |
        | birthday is a 32-bit integer and the test uses 2^12 of them,  |
        | so that j, the number of duplicate spacings, is asympotically |
        | Poisson distributed with lambda=4.  Generators that pass the  |
        | earlier tests for m=1024 and n=2^24 often fail this test, yet |
        | those that pass this test seem to pass the "weaker" test.     |
        | Each set of 4096 birthdays provide a Poisson variate j, and   |
        | 500 such j's lead to a chisquare test to see if the result    |
        | is consistent with the Poisson distribution with lambda=16.   |
        `---------------------------------------------------------------'
            Table of Expected versus Observed counts:
Duplicates    0     1     2     3     4     5     6     7     8     9  >=10
Expected     91   366   732   976   976   781   520   297   148    66    40
Observed     94   395   754  1007   958   740   547   282   126    62    35
(O-E)^2/E   0.1   2.2   0.6   0.9   0.4   2.2   1.3   0.8   3.5   0.3   0.8
            Birthday Spacings: Sum(O-E)^2/E =  13.114, p =  0.784
 ( 4.781 seconds elapsed. )

        .-------------------------------------------------------------.
        | This is the GCD TEST.   Let the (32-bit) RNG produce two    |
        | successive integers u,v.  Use Euclids algorithm to find the |
        | gcd, say x, of u and v. Let k be the number of steps needed |
        | to get x.   Then k is approximately binomial with p=.376    |
        | and n=50,  while the distribution of x is very close to     |
        | Pr(x=i)=c/i^2, with c=6/pi^2.   The gcd test uses ten       |
        | million such pairs u,v to see if the resulting frequencies  |
        | of k's and x's are consistent with the above distributions. |
        | Congruential RNG's---even those with prime modulus---fail   |
        | this test for the distribution of k, the number of steps,   |
        | and often for the distribution of gcd values x as well.     |
        `-------------------------------------------------------------'
Euclid's algorithm:
 p-value, steps to gcd: 0.629783
 p-value, distance of gcd's: 0.126696
 ( 2.086 seconds elapsed. )

        .---------------------------------------------------------------.
        |   This is the GORILLA test, a strong version of the monkey    |
        |   tests that I developed in the 70's. It concerns strings     |
        |   formed from specified bits in 32-bit integers from the RNG. |
        |   We specify the bit position to be studied, from 0 to 31,    |
        |   say bit 3. Then we generate 67,108,889 (2^26+25) numbers    |
        |   from the generator and form a string of 2^26+25 bits by     |
        |   taking bit 3 from each of those numbers. In that string of  |
        |   2^26+25 bits we count the number of 26-bit segments that    |
        |   do not appear. That count should be approximately normal    |
        |   with mean 24687971 and std. deviation 4170. This leads to   |
        |   a normal z-score and hence to a p-value. The test is        |
        |   applied for each bit position 0 (leftmost) to 31.           |
        |   (Some older tests use Fortran's 1-32 for most- to least-    |
        |   significant bits. Gorilla and newer tests use C's 0 to 31.) |
        `---------------------------------------------------------------'
Gorilla test for 2^26 bits, positions 0 to 31:
Note: lengthy test -- ~20 minutes for 900 MHz PC
Bits  0 to  7 || 0.490 0.649 0.276 0.334 0.468 0.222 0.801 0.521
Bits  8 to 15 || 0.619 0.355 0.495 0.452 0.275 0.353 0.053 0.073
Bits 16 to 23 || 0.491 0.109 0.555 0.344 0.552 0.436 0.924 0.169
Bits 24 to 31 || 0.052 0.798 0.589 0.192 0.361 0.438 0.961 0.881
ADKS test for the above 32 p values:  0.751
 ( 27.640 seconds elapsed. )

        .---------------------------------------------------------------.
        |           THE OVERLAPPING 5-PERMUTATION TEST                  |
        | This is the OPERM5 test.  It looks at a sequence of ten       |
        | million 32-bit random integers.  Each set of five consecutive |
        | integers can be in one of 120 states, for the 5! possible     |
        | orderings of five numbers.  Thus the 5th, 6th, 7th,...numbers |
        | each provide a state. As many thousands of state transitions  |
        | are observed,  cumulative counts are made of the number of    |
        | occurences of each state.  Then the quadratic form in the     |
        | weak inverse of the 120x120 covariance matrix yields a test   |
        | that the 120 cellcounts came from the specified (asymptotic)  |
        | distribution with the specified means and 120x120 covariance. |
        `---------------------------------------------------------------'
      The OPERM5 test for 10 million (overlapping) 5-tuples for rand-AvdH.
      p-values for 5 runs: 1.0000 1.0000 1.0000 1.0000 1.0000
 ( 6.984 seconds elapsed. ) ok

[toc] | [prev] | [next] | [standalone]

#26947

Albert van der Horst albert@spenarnc.xs4all.nl wrote:

> Looking into random number generators lately, mhx mentioned the
> diehard program.
> This was written by Marsaglia and with his linux version of the program
> he distributes a file mwcp.ps explaining the Multiply With Carry random
> generators.
> 
> They are simple enough to be typed in by heart whenever you need one,
> yet they have respectable periods of 2^64 for 32-bits and 2^128 for
> 64 bits. Where Marsaglia promotes them they must be pretty good.
> The nice thing are the constants. They are small enough to remember.
> The multipliers are from wikipedia.
> Like the congruentual generators they are simpler in Forth than in other
> languages, especially given the double precision wordsets.
> 
...some [untested] code snipped...

It turns out that both Fortran and C encounter some complexities in 
getting 64-bit randoms from MWC algorithms.  It's certainly simple in 
Forth; you do just what I posted years ago with suitable choices for 
64-bit systems, which I didn't cover way back then.

For possible convenience I'm including a loadable file.  If you really 
want 64-bit randoms on a smaller system, just run separate generators for 
pieces of the value you need.

regards   cgm

\ 64-bit pseudo-random number generators 
\ Charles G. Montgomery   Nov2012  

DEFER rand	\ set to choice of generators

HEX

\  multiply-with-carry (George Marsaglia)
400000000000001 Constant rmult    \ 2^58+1

2VARIABLE rloc
3456123456123456 4567890987654321 rloc 2!   \ Marsaglia's values
: mwc  ( -- u ) rloc 2@ rmult UM* ROT 0 D+ OVER rloc 2! ;

' mwc is rand

\ recommended rmult values for other cell sizes
\ 32 bits: 7A3FFD4B 7549D83B  70DCB7DC  6A364C22  62AFDC7A
\ 16 bits: 65E8 61BF 62DC 6594 6363 5E9B
\  8 bits: 65FF 56FF 68FF 317F
\ periods ~2^29 for 16 bits, ~2^58 for 32 bits, ~2^121 for 64 bits


\  linear congruential multiplier (with Marsaglia's constants)
   
19BAFFBED CONSTANT rgen    12D687 CONSTANT rplus 

VARIABLE rseed   12345678 rseed !
: lcm  ( -- u )  rseed @ rgen um* drop rplus + dup rseed ! ;
\ ====end of file========

[toc] | [prev] | [standalone]

Back to top | Article view | comp.lang.forth

csiph-web