KISS4691, a potentially top-ranked RNG.
in [Fortran]
|
Prev: solutions book
Next: real kind declaration
From: geo on 24 Jul 2010 09:02 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) */
From: Geoff on 24 Jul 2010 14:36 On Sat, 24 Jul 2010 06:02:27 -0700 (PDT), geo <gmarsaglia(a)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) >*/ The compiler also notices that i is defined but unused in MWC(). Purely cosmetic changes for readability and style: /* * 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. * George Marsaglia */ static unsigned long xs = 521288629; static unsigned long xcng = 362436069; static unsigned long Q[4691]; /* takes about 4.2 nanosecs or 238 million/second */ unsigned long MWC(void) { 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); return 0; } /* * 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) */
From: jacob navia on 24 Jul 2010 16:35 geo a �crit : > 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) > */ This doesn't work with systems that have unsigned long as a 64 bit quantity. I obtain: MWC result=3740121002 ? 4169348530 KISS result=2224631993 ? 1421918629 Compiling with 32 bit machine yields: MWC result=3740121002 ? 3740121002 KISS result=2224631993 ? 2224631993
From: Gene on 24 Jul 2010 20:49 On Jul 24, 9:02 am, geo <gmarsag...(a)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;
From: Gib Bogle on 24 Jul 2010 23:34
geo wrote: > > 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) > */ Hi George, I translated this into Fortran, and found that I get different results than with C. I've tracked the difference into MWC(). The following Fortran code, with my laborious comparison of two signed integers treating them as unsigned, gives correct results. If I comment out the line c = tLTx + SHIFTR(x,19) and uncomment the following lines that implement your suggestion above to compute c, I get different results. integer function MWC() integer :: t, x, i integer, save :: c = 0, j = 4691 integer :: tLTx if (j < 4690) then j = j + 1 else j = 0 endif x = Q(j) t = SHIFTL(x,13) + c + x if ((t >= 0 .and. x >= 0) .or. (t < 0 .and. x < 0)) then if (t < x) then tLTx = 1 else tLTx = 0 endif elseif (x < 0) then tLTx = 1 elseif (t < 0) then tLTx = 0 endif c = tLTx + SHIFTR(x,19) !if (sign(1,SHIFTL(x,13)+c) == sign(1,x)) then ! c = sign(1,x) + SHIFTR(x,19) !else ! c = 1 - sign(1,t) + SHIFTR(x,19) !endif Q(j) = t MWC = t end function Is it the case that although your suggested workaround gives different results from the C code in this case, it is still equivalent as a RNG? Cheers Gib |