Re: KISS4691, a potentially top-ranked RNG.

[Gene <gene.ressler@gmail.com>]

On Jul 24, 9:02 am, geo <gmarsag...@gmail.com> wrote:
> I have been asked to recommend an RNG
> (Random Number Generator) that ranks
> at or near the top in all of the categories:
> performance on tests of randomness,
> length of period, simplicity and speed.
> The most important measure, of course, is
> performance on extensive tests of randomness, and for
> those that perform well, selection may well depend
> on those other measures.
>
> The following KISS version, perhaps call it KISS4691,
> seems to rank at the top in all of those categories.
> It is my latest, and perhaps my last, as, at age 86,
> I    am        slowing                down.
>
> Compiling and running the following commented
> C listing should produce, within about four seconds,
> 10^9 calls to the principal component MWC(), then
> 10^9 calls to the KISS combination in another ~7 seconds.
>
> Try it; you may like it.
>
> George Marsaglia
>
> /*
> The KISS4691 RNG, a Keep-It-Simple-Stupid combination of
> MWC (Multiply-With-Carry), Congruential and Xorshift sequences.
> Expressed as a simple MWC() function and C macros
>  #define CNG ( xcng=69069*xcng+123 )
>  #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
>  #define KISS ( MWC()+CNG+XS )
> but easily expressed in other languages, with a slight
> complication for signed integers.
>
> With its immense period, >10^45000, and speed: MWC()s at
> around 238 million/sec or full KISSes at around 138 million,
> there are few RNGs that do as well as this one
> on tests of randomness and are comparable in even one
> of the categories: period, speed, simplicity---not to
> mention comparable in all of them.
>
> The MWC4691 sequence x[n]=8193*x[n-4691]+carry mod b=2^32
> is based on p=8193*b^4691-1, period ~ 2^150124  or 10^45192
> For the MWC (Multiply-With-Carry) process, given a current
> x and c, the new x is    (8193*x+c) mod b,
>          the new c is    the integer part of (8193*x+c)/b.
>
> The routine MWC() produces, in reverse order,  the base-b=2^32
> expansion of some j/p with 0<j<p=8193*2^150112-1 with j
> determined by the 64 bits of seeds xcng and xs, or more
> generally, by 150112 random bits in the Q[] array.
> */
>
> static unsigned long xs=521288629,xcng=362436069,Q[4691];
>
> unsigned long MWC(void)  /*takes about 4.2 nanosecs or 238 million/
> second*/
> {unsigned long t,x,i; static c=0,j=4691;
>   j=(j<4690)? j+1:0;
>   x=Q[j];
>   t=(x<<13)+c+x; c=(t<x)+(x>>19);
>   return (Q[j]=t);
>
> }
>
> #define CNG ( xcng=69069*xcng+123 )
> #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
> #define KISS ( MWC()+CNG+XS ) /*138 million/sec*/
>
> #include <stdio.h>
> int main()
> {unsigned long i,x;
>  for(i=0;i<4691;i++) Q[i]=CNG+XS;
>  for(i=0;i<1000000000;i++) x=MWC();
>  printf(" MWC result=3740121002 ?\n%22u\n",x);
>  for(i=0;i<1000000000;i++) x=KISS;
>  printf("KISS result=2224631993 ?\n%22u\n",x);
>
> }
>
> /*
> This RNG uses two seeds to fill the Q[4691] array by
> means of CNG+XS mod 2^32. Users requiring more than two
> seeds will need to randomly seed Q[] in main().
> By itself, the MWC() routine passes all tests in the
> Diehard Battery of Tests, but it is probably a good
> idea to include it in the KISS combination.
>
> Languages requiring signed integers should use equivalents
> of the same operations, except that the C statement:
>        c=(t<x)+(x>>19);
> can be replaced by that language's version of
>   if sign(x<<13+c)=sign(x) then c=sign(x)+(x>>19)
>            else c=1-sign(t)+(x>>19)
> */

Here's an implementation in Ada, verified to produce the same answers
as the C code.

The Ada version of gcc compiles this to nearly the same code as
produced
for the C. Function call overhead adds about 40% to the run time.

What's a good way to take an arbitrary, single 32-bit seed value to
a complete state initialization?  Accept that lots of people will
pick very small numbers.

-- kiss_random_numbers.ads
-- Specification for George Marsaglia's KISS random number generator.
package KISS_Random_Numbers is

   type Result_Type is mod 2 ** 32;

   type State_Type is private;

   function Next_Random_Value(State : access State_Type) return
Result_Type;

private

   type State_Index_Type is mod 4691;
   type State_Vector_Type is array (State_Index_Type) of Result_Type;
   Initial_Q : State_Vector_Type; -- set in package body

   type Substate_Type is
      record
         Xs   : Result_Type := 521288629;
         Xcng : Result_Type := 362436069;
      end record;

   type State_Type is
      record
         Sub : aliased Substate_Type;
         C   : Result_Type := 0;
         Q   : State_Vector_Type := Initial_Q;
         J   : State_Index_Type := State_Index_Type'Last;
      end record;

end KISS_Random_Numbers;

-- kiss_random_numbers.adb
-- Implementation of George Marsaglia's KISS random number generator.
package body KISS_Random_Numbers is

   function Mwc (State : access State_Type) return Result_Type is
      T, X : Result_Type;
   begin
      State.J := State.J + 1;
      X := State.Q(State.J);
      T := X * 2**13 + State.C + X;
      if T < X then
         State.C := X / 2**19 + 1;
      else
         State.C := X / 2**19;
      end if;
      State.Q(State.J) := T;
      return T;
   end Mwc;
   pragma Inline(Mwc);

   function Cng (State : access Substate_Type) return Result_Type is
   begin
      State.Xcng := 69069 * State.Xcng + 123;
      return State.Xcng;
   end Cng;
   pragma Inline(Cng);

   function Xs (State : access Substate_Type) return Result_Type is
      Xs : Result_Type;
   begin
      Xs := State.Xs;
      Xs := Xs xor (Xs * 2**13);
      Xs := Xs xor (Xs / 2**17);
      Xs := Xs xor (Xs * 2**5);
      State.Xs := Xs;
      return Xs;
   end Xs;
   pragma Inline(Xs);

   function Next_Random_Value(State : access State_Type) return
Result_Type is
   begin
      return Mwc(State) + Cng(State.Sub'Access) +
Xs(State.Sub'Access);
   end Next_Random_Value;

begin
   declare
      S : aliased Substate_Type;
   begin
      for I in Initial_Q'Range loop
         Initial_Q(I) := Cng(S'access) + Xs(S'Access);
      end loop;
   end;
end KISS_Random_Numbers;

-- test_kiss.adb
-- Simple test of George Marsaglia's KISS random number generator.
with Ada.Text_IO; use Ada.Text_IO;
with KISS_Random_Numbers;
use  KISS_Random_Numbers;

procedure Test_Kiss is
   X : Result_Type;
   S : aliased State_Type;
begin
   for I in 1 .. 1000000000 loop
      X := Next_Random_Value(S'Access);
   end loop;
   Put(Result_Type'Image(X));
   New_Line;
end;