Marsaglia's pseudo random number generators

From Lazarus wiki
Revision as of 12:44, 10 March 2017 by Thaddy (talk | contribs)
Jump to navigationJump to search

Marsaglia's pseudo random number generators

Although Freepascal has a reasonably good pseudo random number generator (Prng), a Mersenne Twister, it is rather slow.
Professor George Marsaglia has described a set of Prng's with good properties that are much faster and often just as good.
He was kind enough to provide sourcecode - in C - and a description of each.
Here we present a class that encapsulates some of these for 32 bit random values with a short explanation.
The class contains the following Prng's with Marsaglia's comments from his usenet post. It is extremely fast.
Marsaglia describes them as follows:

multiply with carry

     The MWC generator concatenates two 16-bit multiply-
     with-carry generators, x(n)=36969x(n-1)+carry,
     y(n)=18000y(n-1)+carry mod 2^16, has period about
     2^60 and seems to pass all tests of randomness. A
     favorite stand-alone generator---faster than KISS,
     which contains it.  

fibonaci

     FIB is the classical Fibonacci sequence
     x(n)=x(n-1)+x(n-2),but taken modulo 2^32.
     Its period is 3*2^31 if one of its two seeds is odd
     and not 1 mod 8. It has little worth as a RNG by
     itself, but provides a simple and fast component for
     use in combination generators.
     The classical Fibonacci sequence mod 2^32 from FIB
     fails several tests. It is not suitable for use by
     itself, but is quite suitable for combining with
     other generators.

congruential generator

     CONG is a congruential generator with the widely used 69069
     multiplier: x(n)=69069x(n-1)+1234567. It has period
     2^32. The leading half of its 32 bits seem to pass
     tests, but bits in the last half are too regular.
     The last half of the bits of CONG are too regular,
     and it fails tests for which those bits play a
     significant role. CONG+FIB will also have too much
     regularity in trailing bits, as each does. But keep
     in mind that it is a rare application for which
     the trailing bits play a significant role. CONG
     is one of the most widely used generators of the
     last 30 years, as it was the system generator for
     VAX and was incorporated in several popular
     software packages, all seemingly without complaint.
     (note: this is about the same as Delphi uses)

four lagged Fibonaci

     LFIB4 is an extension of what (Marsaglia) have previously
     defined as a lagged Fibonacci generator:
     x(n)=x(n-r) op x(n-s), with the x's in a finite
     set over which there is a binary operation op, such
     as +,- on integers mod 2^32, * on odd such integers,
     exclusive-or(xor) on binary vectors. Except for
     those using multiplication, lagged Fibonacci
     generators fail various tests of randomness, unless
     the lags are very long. (See SWB below).
     To see if more than two lags would serve to overcome
     the problems of 2-lag generators using +,- or xor, I
     have developed the 4-lag generator LFIB4 using
     addition: x(n)=x(n-256)+x(n-179)+x(n-119)+x(n-55)
     mod 2^32. Its period is 2^31*(2^256-1), about 2^287,
     and it seems to pass all tests---in particular,
     those of the kind for which 2-lag generators using
     +,-,xor seem to fail. For even more confidence in
     its suitability, LFIB4 can be combined with KISS,
     with a resulting period of about 2^410: just use
     (KISS+LFIB4) in any Pascal expression.

3-shift-register generator

     SHR3 is a 3-shift-register generator with period
     2^32-1. It uses y(n)=y(n-1)(I+L^17)(I+R^13)(I+L^5),
     with the y's viewed as binary vectors, L the 32x32
     binary matrix that shifts a vector left 1, and R its
     transpose. SHR3 seems to pass all except those
     related to the binary rank test, since 32 successive
     values, as binary vectors, must be linearly
     independent, while 32 successive truly random 32-bit
     integers, viewed as binary vectors, will be linearly
     independent only about 29% of the time.
     (note: this is also known as XORSHIFT and very fast)

shift with borrow

     SWB is a subtract-with-borrow generator that I
     developed to give a simple method for producing
     extremely long periods:
     x(n)=x(n-222)-x(n-237)- borrow mod 2^32.
     The 'borrow' is 0, or set to 1 if computing x(n-1)
     caused overflow in 32-bit integer arithmetic. This
     generator has a very long period, 2^7098(2^480-1),
     about 2^7578. It seems to pass all tests of
     randomness, except for the Birthday Spacings test,
     which it fails badly, as do all lagged Fibonacci
     generators using +,- or xor. I would suggest
     combining SWB with KISS, MWC, SHR3, or CONG.
     KISS+SWB has period >2^7700 and is highly
     recommended.
     Subtract-with-borrow has the same local behaviour
     as lagged Fibonacci using +,-,xor---the borrow
     merely provides a much longer period.
     SWB fails the birthday spacings test, as do all
     lagged Fibonacci and other generators that merely
     combine two previous values by means of =,- or xor.
     Those failures are for a particular case: m=512
     birthdays in a year of n=2^24 days. There are
     choices of m and n for which lags >1000 will also
     fail the test. A reasonable precaution is to always
     combine a 2-lag Fibonacci or SWB generator with
     another kind of generator, unless the generator uses
     *, for which a very satisfactory sequence of odd
     32-bit integers results.

keep it simple sir

     The KISS generator, (Keep It Simple Stupid), is
     designed to combine the two multiply-with-carry
     generators in MWC with the 3-shift register SHR3 and
     the congruential generator CONG, using addition and
     exclusive-or. Period about 2^123.
     It is one of (Marsaglia's) favorite generators.

implementation

The class is a fully static sealed class so you do not have to instantiate it or derive from it.

unit uMarsaglia;
{$mode objfpc}
interface

type
  { It is not meant to be extended, since it is based on a single article so I sealed the class }
  TMarsagliaPrngs = class sealed 
  strict private
    {static variables:}
    class var T: array [0..255] of dword;
    class var a,b,v,w,x,y,z,bro,jcong,jsr:dword;
    class var c:byte;
    class function wnew:dword;inline;static;    
    class function znew:dword;inline;static;
  public
    { You can also call the create explicitly to reset the static 
      variables.}
    class constructor create;
    class procedure settable(const i1,i2,i3,i4,i5,i6:dword);static;
    class procedure Randomize;inline;static;
    class function MWC:dword;inline;static;  
    class function SHR3:dword;inline;static;  
    class function CONG:dword;inline;static;  
    class function FIB:dword;inline;static;  
    class function KISS:dword;inline;static; 
    class function LFIB4:dword;inline;static;
    class function SWB:dword;inline;static;  
  end;
  
  { a shorthand alias: }
  Prng = type TMarsagliaPrngs;     
  
  { Because many simulations call for uniform
    random variables in 0<x<1 or -1<x<1, here's a utility function.
    Using Signed = false will provide a random single in range (0,1), 
    while Signed true will provide one in range (-1,1).
    Simply feed it the output of one of the Prng's}
  function ScaleToFloat(const value:dword; Signed:Boolean=true):single;inline;
     

implementation
  
    class function TMarsagliaPrngs.znew:dword;inline;static;
    begin
      z:=36969*(z and 65535)+(z >> 16);
      Result:= z;
    end;
        
    class function TMarsagliaPrngs.wnew:dword;inline;static;
    begin
      w:=18000*(w and 65535)+(w >> 16);
      Result := w;
    end;
    
    
    class constructor TMarsagliaPrngs.create;
    begin
      a:=224466889; b:=7584631;c:=0;v:=2463534242;w:=521288629;x:=0;y:=0;
      z:=362436069;jsr:=123456789; jcong:=380116160;  
      settable(12345,65435,34221,12345,9983651,95746118);  
    end;
    
    { Use random seeds to reset z,w,jsr,jcong,a,b, and the table t[256]}  
    class procedure TMarsagliaPrngs.settable(const i1,i2,i3,i4,i5,i6:dword);inline;static;
    var
      i:integer;
    begin
      z:=i1; w:=i2; jsr:=i3; jcong:=i4; a:=i5; b:=i6;
    { Example procedure to set the table using KISS:}
      for i:=0 to 255 do t[i]:=KISS;
    end;
    
    class function TMarsagliaPrngs.MWC:dword;inline;static;
    begin
      Result:= (znew shl 16)+wnew; 
    end;
    
    class function TMarsagliaPrngs.SHR3:dword;inline;static;
    begin
      jsr:=jsr xor(jsr << 17); jsr:= jsr xor(jsr >> 13); jsr:=jsr xor(jsr << 5);
      result := jsr;
    end;
    
    class function TMarsagliaPrngs.CONG:dword;inline;static;
    begin
     jcong:=69069*jcong+1234567;
     Result := jcong;
    end;
    
    class function TMarsagliaPrngs.FIB:dword;inline;static;
    begin
       b:=a+b;a:=b-a;
       Result := a;
    end;
    
    class function TMarsagliaPrngs.KISS:dword;inline;static;
    begin
      Result:= (MWC xor CONG)+SHR3;
    end;
    
    class function TMarsagliaPrngs.LFIB4:dword;inline;static;
    begin
      inc(c);
      t[c]:=t[c]+t[(c+58) and 255]+t[(c+119) and 255]+t[(c+178) and 255];
      Result:=t[c];
    end;
    
    class function TMarsagliaPrngs.SWB:dword;inline;static;
    begin
      inc(c);bro:=dword(x<y);
      x:=t[(c+34) and 255];
      y:=t[(c+19) and 255]+bro;
      t[c] := x-y;
      Result:= t[c];
    end;
    
    class procedure TMarsagliaPrngs.Randomize;inline;static;
    begin
      {$WARNING This is not implemented here yet}
    end;
    
   function ScaleToFloat(const value:dword; Signed:Boolean=true):single;inline;
   begin
     case Signed of
     False:Result:= dword(value) * 2.328306e-10;
     True :Result:= longint(value) * 4.656613e-10;
     end;
  end;
     
 end.

example

The original C code contained an example that is here translated to Object Pascal. It is just a proof that we get the same results as intended.

program marsagliatest;
{$mode objfpc}
{ This is a test main program. It should compile and print 8 TRUE's. 
  That means it renders the same values as Marsaglia intended, given the
  seeds provided. It is basically Marsaglia's C test program in Pascal }
uses 
  uMarsaglia;
var
  i:integer;
  k:dword;
begin 
  // set the table to known values, so we can verify the results.
  // note the extreme lag for the FPC random.... 
  Prng.settable(12345,65435,34221,12345,9983651,95746118);
  for i:=1 to 1000000 do k:=Prng.LFIB4;writeln('LFIB4 ','Pass: ',k-1064612766=0);
  for i:=1 to 1000000 do k:=Prng.SWB;  writeln('SWB   ','Pass: ',k- 627749721=0);
  for i:=1 to 1000000 do k:=Prng.KISS; writeln('KISS  ','Pass: ',k-1372460312=0);
  for i:=1 to 1000000 do k:=Prng.CONG; writeln('CONG  ','Pass: ',k-1529210297=0);
  for i:=1 to 1000000 do k:=Prng.SHR3; writeln('SHR3  ','Pass: ',k-2642725982=0);
  for i:=1 to 1000000 do k:=Prng.MWC;  writeln('MWC   ','Pass: ',k- 904977562=0);
  for i:=1 to 1000000 do k:=Prng.FIB;  writeln('FIB   ','Pass: ',k-3519793928=0);
  for i:=1 to 1000000 do k:=Random(high(Dword));writeln('FPC   ','Pass: ',k-1314949284=0);
  // Scale to float demo 
  for i := 0 to 9 do writeln(ScaleToFloat(Prng.SHR3):2:10);
  writeln('..........................');
  for i := 0 to 9 do writeln(ScaleToFloat(Prng.SHR3, false):2:10);
end.