> /* > 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;
> 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 */
/* Here is a C++ version. The C version is quite a bit faster because there are no function calls at all. Can any of you C++ gurus bump the speed without losing encapsulation? I get about 5 seconds for the C version and about 8 seconds for the C++ version.
-- d.corbit */
#include <iostream> /* 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: */
class SuperKiss {
private: unsigned long Q[41790]; unsigned long indx; unsigned long carry; unsigned long xcng; unsigned long 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 int x; SuperKiss sk; for (i = 0; i < 1000000000; i++) x = sk.SKRand(); std::cout << " x = " << x << std::endl << "does Does x=-872412446?" << std::endl; return 0;
}
/* 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
int main () { unsigned long i; int x=0; SuperKiss sk; for (i = 0; i < 1000000000; i++) x = sk.SKRand(); std::cout << " x = " << x << std::endl << "Does x=-872412446?" << std::endl; return 0;
On Nov 4, 1:26 pm, user923005 <dcor...@connx.com> wrote:
> On Nov 4, 9:18 am, David <david.astgt...@gmail.com> wrote:
> > On an x86-64 machine using GCC version 4.3.3 (Ubuntu 4.3.3-5ubuntu4), > > both the C code and C++ code fail for me. > > I get: > > x=505478909. > > Does x=-872412446?
> > Changing the unsigned long's to unsigned int's fixed the problem. > > And it does matter: before the change, the generator failed a variety > > of tests (really odd assortment, though: parking lot, 2dsphere, > > 3dsphere, squeeze, and sums).
> OK, makes sense. The RNG must assume 32 bit longs.
Modified C++ code:
#include <iostream>
class SuperKiss {
private: unsigned int Q[41790]; unsigned int indx; unsigned int carry; unsigned int xcng; unsigned int 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] =(unsigned int) ~(t); } indx = 1; return (Q[0]); }
int main () { unsigned int i; int x=0; SuperKiss sk; for (i = 0; i < 1000000000; i++) x = sk.SKRand(); std::cout << " x = " << x << std::endl << "Does x=-872412446?" << std::endl; return 0;
}
/* Possible output: x = -872412446 Does x=-872412446? */
>I copied and pasted from the wrong file. Here is the correct code >[snip] >class SuperKiss {
this how i see it 1: but it is not much portable (cpu x86, assembler nasm, c++ compiler borland) then don't know if it is right what about this? --xx Delta=20.000 i=-872412446 Does i=-872412446? val u32 =279941939 val double =0.600837 val u32 =1357108688 val i32 -10309 +309 =-6962 val d -10309 +309 =-695.931958 --xx #include <stdio.h> #include <stdint.h> #include <time.h>
#define u32 uint32_t #define i32 int32_t #define R return #define P printf #define F for
class SuperKiss { public: SuperKiss(); /* def constructor */ u32 SKRand(void); /* def function */ u32 urnd(void){R SKRand(); } i32 irnd(void){R (i32) SKRand(); } u32 urnd(u32 umin, u32 umax) {if(umin>=umax) R 0; R umin+SKRand()%(umax-umin+1); } i32 irnd(i32 imin, i32 imax) {if(imin>=imax) R 0; R imin+SKRand()%(imax-imin+1); } double drnd(void){R (double)SKRand()/0xFFFFFFFF;} double drnd(double dmin, double dmax) {if(dmin>=dmax) R 0.0; R dmin + (dmax-dmin)*drnd(); }
};
int main(void) {double d; time_t ti,tf; SuperKiss sk, w; i32 i; u32 u, v;
>I copied and pasted from the wrong file. Here is the correct code >[snip] >class SuperKiss {
this is how i see it 2: but it is not much portable (cpu x86, assembler nasm, c++ compiler borland) then don't know if it is right... what about this? ----x Delta=23.000 i=-872412446 Does i=-872412446? val u32 =3075790285 val double =0.647830 val u32 =4289339058 val i32 -10309 +309 =-6324 val d -10309 +309 =-6716.981099 ---x #include <stdio.h> #include <stdint.h> #include <time.h> #include <stdlib.h>
#define u32 uint32_t #define i32 int32_t #define R return #define P printf #define F for
class SuperKiss{ public: SuperKiss(u32, u32, u32); /* def costructor */ ~SuperKiss(){ free(q); } u32 SKRand(void); /* def function */ u32 urnd(void){R SKRand(); } i32 irnd(void){R (i32) SKRand(); } u32 urnd(u32 umin, u32 umax) {if(umin>=umax) R 0; R umin+SKRand()%(umax-umin+1); } i32 irnd(i32 imin, i32 imax) {if(imin>=imax) R 0; R imin+SKRand()%(imax-imin+1); } double drnd(void){R (double)SKRand()/0xFFFFFFFF;} double drnd(double dmin, double dmax) {if(dmin>=dmax) R 0.0; R dmin + (dmax-dmin)*drnd(); } u32* q; u32 indx; u32 carry; u32 xcng; u32 xs;
|
