RNGs: A Super KISS
geo
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 */
user923005
/*
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>
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]);
}
public:
// Constructor:
SuperKiss()
{
indx = 41790;
carry = 362436;
xcng = 1236789;
xs = 521288629;
unsigned i;
for (i = 0; i < 41790; i++)
Q[i] = (xcng = 69609 * xcng + 123) +
(xs ^= xs << 13, xs ^= (unsigned) xs >> 17, xs ^=
xs >> 5);
}
// Collect next random number:
unsigned long SKRand() {
return (indx < 41790 ? Q[indx++] : refill ()) +
(xcng = 69609 * xcng + 123) +
(xs ^= xs << 13, xs ^= (unsigned) xs >> 17, xs ^= xs >>
5);
}
};
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;
Tom St Denis
> 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]);
> }
Not to nitpick but your C code could use some work. First off, some
indentation please? Second, returning a unsigned long long as "int"
is not very portable.
Part of the good thing of PRNGs is that they're reproduceable.
Ideally over different platforms. You should truncate the return
value if you want it as "int" or change the return type. As it stands
now this will produce different results on my x86-32 and 64 boxes for
the same seed, which is a bad thing.
Tom
user923005
I copied and pasted from the wrong file. Here is the correct code
[snip]
class SuperKiss {
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;
}
James Waldby
> 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.
> static unsigned long Q[41790],indx=41790,carry=362436,xcng=1236789,xs=521288629;
> {unsigned long i,x;
...
> printf(" x=%d.\nDoes x=-872412446?\n",x);
> Running this program should produce 10^9 KISSes in some 7-15 seconds.
> 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
I've snipped the program except for three lines that apparently must
differ depending upon cpu word length. On my 64-bit Athlon X2 5200+
(1GHz) with gcc 4.1.2, and %d changed to %ld, the output (after 7.5
seconds) contains "x=2904265093743181565."; or, with instead long
changed to int in two places, "x=-872412446." (after 7.3 seconds).
--
jiw
James Waldby
Note, in hexadecimal those results are 284e05b71e20fefd and cc000ae2
respectively. Ie, the bit patterns of the low 32 bits are quite
different.
--
jiw
user923005
I get the same results on:
64 bit Windows using the 64 bit MS compiler
64 bit Windows using the 32 bit MS compiler
64 bit Windows using the 64 bit Mingw GCC compiler
64 bit OpenVMS (Itanium) using HP CXX
64 bit OpenVMS (Alpha) using HP CXX
32 bit OpenVMS (VAX) using HP CXX (Had to remove the std:: because of
old compiler, expect a very long wait)
Solaris 5.9 is interesting because it is big-endian, in contrast with
those previously mentioned:
/export/home/dcorbit> uname -a
SunOS solaris9 5.9 Generic_118558-11 sun4u sparc SUNW,Sun-Fire-V210
/export/home/dcorbit> gcc --version
gcc (GCC) 4.0.2
Copyright (C) 2005 Free Software Foundation, Inc.
This is free software; see the source for copying conditions. There
is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE.
/export/home/dcorbit> ./a.out
x = -872412446
Does x=-872412446?
/export/home/dcorbit>
Dann Corbit
@d21g2000yqn.googlegroups.com>, t...@iahu.ca says...
>
> On Nov 3, 10:46ï¿œam, geo <gmarsag...@gmail.com> wrote:
> > ￜ 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]);
> > ￜ }
>
> Not to nitpick but your C code could use some work. First off, some
> indentation please? Second, returning a unsigned long long as "int"
> is not very portable.
The O.P. is George Marsaglia. If there were a Mt. Rushmore for computer
science, Donald Knuth would be George Washington's bust, but Marsaglia
would be up there somewhere too. If it were random numbers, then Mr.
Marsaglia is front and center.
IMO-YMMV
> Part of the good thing of PRNGs is that they're reproduceable.
> Ideally over different platforms. You should truncate the return
> value if you want it as "int" or change the return type. As it stands
> now this will produce different results on my x86-32 and 64 boxes for
> the same seed, which is a bad thing.
Did you actually try it?
What compilers were you using?
I get the same result regardless of compiler, hardware and OS.
io_x
"geo" <gmars...@gmail.com> ha scritto nel messaggio
news:a14dbdd8-cce7-4038...@f16g2000yqm.googlegroups.com...
> #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);
> }
the asm version
; compile with
; nasmw -fobj rndasm.asm
; bcc32 -v rndasm.obj
section _DATA use32 public class=DATA
global _main
extern _printf
indx dd 41790
carry dd 362436
xcng dd 1236789
xs dd 521288629
IIIIIIxIIdIIn db " x=%d." , 13, 10, 0
IDoesIxII872412446IIn db "Does x=-872412446?" , 13, 10, 0
section _BSS use32 public class=BSS
vettoreQ resd 41792
section _TEXT use32 public class=CODE
align 4
initialize:
push esi
push edi
push ebp
mov ecx, 41790
mov esi, vettoreQ
mov edi, [xcng]
mov ebp, [xs]
.0: mov eax, 69609
mul edi
add eax, 123
xchg eax, edi
mov eax, ebp
shl eax, 13
xor ebp, eax
mov eax, ebp
shr eax, 17
xor ebp, eax
mov eax, ebp
shr eax, 5
xor ebp, eax
lea eax, [ebp+edi]
mov [esi], eax
add esi, 4
loop .0
mov [xcng], edi
mov [xs], ebp
pop ebp
pop edi
pop esi
ret
align 4
rnd:
push esi
push edi
mov ecx, [indx]
mov esi, vettoreQ
cmp ecx, 41790
jne .1
push esi
.0: mov eax, 7010176
mul dword[esi]
add eax, [carry]
adc edx, 0
mov [carry], edx
not eax
mov dword[esi], eax
add esi, 4
loop .0
pop esi
.1: mov edi, [esi+4*ecx]
inc ecx
mov [indx], ecx
mov eax, 69609
mov ecx, [xcng]
mul ecx
add eax, 123
mov [xcng], eax
add edi, eax
mov edx, [xs]
mov eax, edx
shl eax, 13
xor edx, eax
mov eax, edx
shr eax, 17
xor edx, eax
mov eax, edx
shr eax, 5
xor edx, eax
mov [xs], edx
add edi, edx
xchg edi, eax
pop edi
pop esi
ret
align 4
_main:
pushad
call initialize
mov esi, 0
.0: call rnd
inc esi
cmp esi, 1000000000
jb .0
mov edi, eax
push eax
push IIIIIIxIIdIIn
call _printf
add esp, 8
push IDoesIxII872412446IIn
call _printf
add esp, 4
popad
xor eax, eax
ret
io_x
"Tom St Denis" <t...@iahu.ca> ha scritto nel messaggio
news:4fd3c8bf-e45b-4359...@d21g2000yqn.googlegroups.com...
On Nov 3, 10:46 am, geo <gmarsag...@gmail.com> wrote:
> 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]);
> }
Not to nitpick but your C code could use some work. First off, some
indentation please? Second, returning a unsigned long long as "int"
is not very portable.
<where is written that he return long long int like int??
<Q is "static unsigned long Q[41790]"
<he return only "unsigned long" like "int"
<the error could be here "Q[i]=~(t)" because in the left side is long
<the other side is long long
David
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).
David
James Waldby
> 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?
That result is like what I mentioned in my posts yesterday;
505478909. == 0x1e20fefd, which is the low 32 bits
of 2904265093743181565. == 0x284e05b71e20fefd .
> 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).
...
--
jiw
user923005
OK, makes sense. The RNG must assume 32 bit longs.
user923005
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]);
}
public:
// Constructor:
SuperKiss()
{
indx = 41790;
carry = 362436;
xcng = 1236789;
xs = 521288629;
unsigned i;
for (i = 0; i < 41790; i++)
Q[i] = (xcng = 69609 * xcng + 123) +
(xs ^= xs << 13, xs ^= (unsigned) xs >> 17, xs ^=
xs >> 5);
}
// Collect next random number:
unsigned int SKRand() {
return (indx < 41790 ? Q[indx++] : refill ()) +
(xcng = 69609 * xcng + 123) +
(xs ^= xs << 13, xs ^= (unsigned) xs >> 17, xs ^= xs >>
5);
}
};
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?
*/
Keith
> 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.
I've taken a look through this, and it's not clear to me how you 'seed'
this generator to provide a unique 'starting point'. Is it sufficient
to initiallise Q[] with random stuff?
geo
Yes.
As stated in the original posting:
"Users should choose a reasonable number of the 1337280
random bits that a fully-seeded Q array requires."
The degree of reasonableness depends on whether the
application is for simulation or for security.
With b=2^32, the SUPR part of this KISS=SUPR+CNG+XS
produces, in reverse order by the CMWC method, the
base-b expansion of the rational number k/p for
the prime p=7010176*b^41490+1, with k determined
by the seed values in Q and the initial carry.
There are repeating cycles of 54767*2^1337279 base-b
"digits" in that expansion---the order of b for the prime p.
If you have a given set of a few thousand "digits" in
that sequence and you want to know precursors or
successors, you are pretty much stuck if the initial
Q array was fully seeded---there are so many places
where a particular string can appear that finding a
specific one becomes extremely difficult.
If only, say, 64 bits were used to seed the Q array
by means of CNG+XS, then a particular string may appear
at considerably fewer places in that base-b expansion.
The task of locating a particular location, so as to
get precursors or successors, is reduced but may be feasible.
But any one of the many appearances of that particular
string of base-b digits is likely to serve well as a set
of independent random 32-bit integers---just, as, for example,
strings of bits in the binary expansion of pi are could well
run Las Vegas, if chosen from an unknown random location
among the first 10^400000 places, specified each day.
Thus, for simulation purposes, some 60-100 seed bits seem to
serve very well; for security, more of the 1.3 million
seed bits should chosen---the more, the safer.