Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: RNGs: A Super KISS
Replies: 14   Last Post: Feb 26, 2013 10:53 AM

 Messages: [ Previous | Next ]
 David Bernier Posts: 3,892 Registered: 12/13/04
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 Mersenne-ized?) by a RNG
>>> with period 2^19937-1, I offer one here with period
>>> 54767*2^1337279---over 10^396564 times as long.
>>> It is one of my CMWC (Complimentary-Multiply-With-Carry) RNGs,
>>> and is suggested here as one of the components of a
>>> super-long-period KISS (Keep-It-Simple-Stupid) RNG.
>>>
>>> With b=2^32 and a=7010176, and given a 32-bit x, and a 32-bit c, this
>>> generator produces a new x,c by forming 64-bit t=a*x+c then replacing:
>>> c=top 32 bits of t and x=(b-1)-(bottom 32 bits of t). In C: c=t>>32;
>>> x=~t;
>>>
>>> For many years, CPUs have had machine instructions to form such a
>>> 64-bit 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 fast---some 140 million per second on my desktop PC.
>>>
>>> Used alone, this generator passes all the Diehard Battery of Tests,
>>> but
>>> its simplicity makes it well-suited to serve as one of the three
>>> components
>>> of a KISS RNG, based on the Keep-It-Simple-Stupid principle, and the
>>> idea,
>>> supported by both theory and practice, that the combination of RNGs
>>> based on
>>> different mathematical models can be no worse---and is usually
>>> better---than
>>> 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 super-long-period 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 7-15 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 (64-bit)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 base-b 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 base-b 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 fully-seeded 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

Marsaglia wrote:
m=(2^28-1)b^(2^22)-1=(2^28-1)b^4194304-1 with b=2^32,
m=(2^28-1)B^(2^21)-1=(2^28-1)B^2091752-1 with B=2^64,

with b=2^32, B=2^64.

2^21 = 2097152 so in fact,

m=(2^28-1)b^(2^22)-1=(2^28-1)b^4194304-1 with b=2^32,
m=(2^28-1)B^(2^21)-1=(2^28-1)B^2097152-1 with B=2^64,

(same m).

In 64-bit, or base 2^64, the state vector has about
2,097,152 64-bit numbers.

MWC is "multiply with carry".

"when put in reverse order, will be the base-b expansion of a rational
j/(ab^r ? 1) for some 0 < j < ab^r."

cf.:
< http://en.wikipedia.org/wiki/Multiply-with-carry >

r is very large, so the lag is large.

Anyway, it works at 67 million 64-bit 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.

Date Subject Author
11/3/09 geo
11/3/09 Dann Corbit
11/3/09 tom@iahu.ca
11/3/09 Dann Corbit