
Re: RNGs: A Super KISS
Posted:
Feb 26, 2013 10:53 AM


On 02/26/2013 08:06 AM, rossum wrote: > On Mon, 25 Feb 2013 20:10:54 0800 (PST), latexpst@gmail.com wrote: > >> On Wednesday, November 4, 2009 2:46:27 AM UTC+11, geo wrote: >>> /* >>> For those mesmerized (or Mersenneized?) by a RNG >>> with period 2^199371, I offer one here with period >>> 54767*2^1337279over 10^396564 times as long. >>> It is one of my CMWC (ComplimentaryMultiplyWithCarry) RNGs, >>> and is suggested here as one of the components of a >>> superlongperiod KISS (KeepItSimpleStupid) RNG. >>> >>> With b=2^32 and a=7010176, and given a 32bit x, and a 32bit c, this >>> generator produces a new x,c by forming 64bit t=a*x+c then replacing: >>> c=top 32 bits of t and x=(b1)(bottom 32 bits of t). In C: c=t>>32; >>> x=~t; >>> >>> For many years, CPUs have had machine instructions to form such a >>> 64bit t and extract the top and bottom halves, but unfortunately >>> only recent Fortran versions have means to easily invoke them. >>> >>> Ability to do those extractions leads to implementations that are >>> simple >>> and extremely fastsome 140 million per second on my desktop PC. >>> >>> Used alone, this generator passes all the Diehard Battery of Tests, >>> but >>> its simplicity makes it wellsuited to serve as one of the three >>> components >>> of a KISS RNG, based on the KeepItSimpleStupid principle, and the >>> idea, >>> supported by both theory and practice, that the combination of RNGs >>> based on >>> different mathematical models can be no worseand is usually >>> betterthan >>> any of the components. >>> >>> So here is a complete C version of what might be called a SUPER KISS >>> RNG, >>> combining, by addition mod 2^32, a Congruential RNG, a Xorshift RNG >>> and the superlongperiod CMWC RNG: >>> */ >>> >>> #include<stdio.h> >>> static unsigned long Q >>> [41790],indx=41790,carry=362436,xcng=1236789,xs=521288629; >>> >>> #define CNG ( xcng=69609*xcng+123 ) /*Congruential*/ >>> #define XS ( xs^=xs<<13, xs^=(unsigned)xs>>17, xs^=xs>>5 ) / >>> *Xorshift*/ >>> #define SUPR ( indx<41790 ? Q[indx++] : refill() ) >>> #define KISS SUPR+CNG+XS >>> >>> int refill( ) >>> { int i; unsigned long long t; >>> for(i=0;i<41790;i++) { t=7010176LL*Q[i]+carry; carry=(t>>32); Q[i]=~ >>> (t);} >>> indx=1; return (Q[0]); >>> } >>> >>> int main() >>> {unsigned long i,x; >>> for(i=0;i<41790;i++) Q[i]=CNG+XS; >>> for(i=0;i<1000000000;i++) x=KISS; >>> printf(" x=%d.\nDoes x=872412446?\n",x); >>> } >>> >>> /* >>> Running this program should produce 10^9 KISSes in some 715 seconds. >>> You are invited to cut, paste, compile and run for yourself, checking >>> to >>> see if the last value is as designated, (formatted as a signed integer >>> for >>> potential comparisons with systems using signed integers). >>> You may want to report or comment on implementations for other >>> languages. >>> >>> The arithmetic operations are suited for either signed or unsigned >>> integers. >>> Thus, with (64bit)t=a*x+c, x=t%b in C or x=mod(t,b) in Fortran, and >>> c=c/b in either C or Fortran, but with ways to avoid integer >>> divisions, >>> and subsequent replacement of x by its baseb complement, ~x in C. >>> >>> With b=2^32 and p=54767*2^1337287+1, the SUPR part of this Super KISS >>> uses my CMWC method to produce, in reverse order, the baseb expansion >>> of k/p for some k determined by the values used to seed the Q array. >>> The period is the order of b for that prime p: >>> 54767*2^1337279, about 2^1337294 or 10^402566. >>> (It took a continuous run of 24+ days on an earlier PC to >>> establish that order. My thanks to the wizards behind PFGW >>> and to Phil Carmody for some suggested code.) >>> >>> Even the Q's all zero, should seeding be overlooked in main(), >>> will still produce a sequence of the required period, but will >>> put the user in a strange and exceedingly rare place in the entire >>> sequence. Users should choose a reasonable number of the 1337280 >>> random bits that a fullyseeded Q array requires. >>> >>> Using your own choices of merely 87 seed bits, 32 each for xcng,xs >>> and 23 for carry<7010176, then initializing the Q array with >>> for(i=0;i<41790;i++) Q[i]=CNG+XS; >>> should serve well for many applications, but others, such as in >>> Law or Gaming, where a minimum number of possible outcomes may be >>> required, might need more of the 1337280 seed bits for the Q array. >>> >>> As might applications in cryptography: With an unknown but fully >>> seeded Q array, a particular string of, say, 41000 successive SUPR >>> values will appear at more than 2^20000 locations in the full >>> sequence, >>> making it virtually impossible to get the location of that particular >>> string in the full loop, and thus predict coming or earlier values, >>> even if able to undo the CNG+XS operations. >>> */ >>> >>> /* >>> So I again invite you to cut, paste, compile and run the above C >>> program. >>> 1000 million KISSes should be generated, and the specified result >>> appear, >>> by the time you count slowly to fifteen. >>> (Without an optimizing compiler, you may have to count more slowly.) >>> */ >>> >>> /* George Marsaglia */ >> >> What's the best way to generate double precision numbers between 0 and 1, excluding 0, from this generator? Thanks. > You may have to wait some time for a reply. Unfortunately George died > in 2011: http://en.wikipedia.org/wiki/George_Marsaglia > > rossum
There was another post by Marsaglia in January 2011:
http://mathforum.org/kb/plaintext.jspa?messageID=7359611
and treeview of thread: http://mathforum.org/kb/thread.jspa?threadID=2228415&messageID=7359611#replytree
Marsaglia wrote: m=(2^281)b^(2^22)1=(2^281)b^41943041 with b=2^32, m=(2^281)B^(2^21)1=(2^281)B^20917521 with B=2^64,
with b=2^32, B=2^64.
2^21 = 2097152 so in fact,
m=(2^281)b^(2^22)1=(2^281)b^41943041 with b=2^32, m=(2^281)B^(2^21)1=(2^281)B^20971521 with B=2^64,
(same m).
In 64bit, or base 2^64, the state vector has about 2,097,152 64bit numbers.
MWC is "multiply with carry".
"when put in reverse order, will be the baseb expansion of a rational j/(ab^r ? 1) for some 0 < j < ab^r."
cf.: < http://en.wikipedia.org/wiki/Multiplywithcarry >
r is very large, so the lag is large.
Anyway, it works at 67 million 64bit unsigned long longs per second:
[david@localhost next]$ gcc testnew01a.c
[david@localhost next]$ ./a.out Does x=13596816608992115578 ? x=13596816608992115578 Does x=5033346742750153761 ? x=5033346742750153761
[david@localhost next]$ time ./a.out Does x=13596816608992115578 ? x=13596816608992115578 Does x=5033346742750153761 ? x=5033346742750153761
real 0m29.444s user 0m29.335s sys 0m0.014s
David Bernier
 dracut:/# lvm vgcfgrestore File descriptor 9 (/.console_lock) leaked on lvm invocation. Parent PID 993: sh Please specify a *single* volume group to restore.

