Re: SuperKISS for 32- and 64-bit RNGs in PL/I
robin
| On Nov 3 I posted
|
| RNGs: A Super KISS sci.math, comp.lang.c, sci.crypt
|
| a KISS (Keep-It-Simple-Stupid) RNG combining,
| by addition mod 2^32, three simple RNGs:
| CMWC(Complementary-Multiply-With-Carry)+CNG(Congruential)
| +XS(Xorshift)
| with resulting period greater than 10^402575.
|
| The extreme period comes from finding a prime p=a*b^r+1 for
| which the order of b has magnitude near p-1, then use
| the CMWC method, the mathematics of which I outline here:
|
| Let Z be the set of all "vectors" of the form
| [x_1,...,x_r;c] with 0<=x_i<b and 0<=c<a.
| Then Z has ab^r elements, and if the function f() on Z is
| f([x_1,x2,...,x_r;c])=
| [x_2,...,x_r,b-1-(t mod b);trunc((t/b)], t=a*x_1+c
| then f() has an inverse on Z: for each z in Z there is exactly
| one w in Z for which f(w)=z:
| f^{-1}([x_1,x2,...,x_r;c])=
| [trunc((v/a),x_1,...,x_{r-1}; v mod a], v=cb+(b-1)-x_r
|
| Thus a directed graph based on z->f(z) will consist only
| of disjoint loops of size s=order(b,p) and there will be
| L=ab^r/s such loops.
|
| A random uniform choice of z from Z is equally likely to
| fall in any one of the L loops, and then each "vector" in
| the sequence obtained by iterating f:
| f(z), f(f(z)), f(f(f(z))),...
| will have a uniform distribution over that loop, and the sequence
| determined by taking the r'th element from each "vector",
| (the output of the CMWC RNG) will be periodic with period
| the order of b for the prime p=ab^r+1, and that sequence
| will produce, in reverse order, the base-b expansion of a
| rational k/p for some k determined by choice of the seed z.
|
| For that Nov 3 post, I had found that the order of b=2^32
| for the prime p=7010176*b^41790+1 is 54767*2^1337279, about
| 10^402566, thus providing an easily-implemented
| KISS=CMWC+CNG+XS RNG with immense period and
| excellent performance on tests of randomness.
|
| That easy implementation required carrying out the essential
| parts of the CMWC function f(): form t=7010176*x+c in 64 bits,
| extract the top 32 bits for the new c, the bottom 32 for
| the new x---easy to do in C, not easy in Fortran.
| And if we want 64-bit random numbers, with B=2^64, our prime
| becomes 7010176*B^20985+1, for which the period of B is
| 54767*2^1337278, still immense, but in C, with 64-bit x,c
| there seems no easy way to form t=7010176*x+c in 128 bits,
| then extract the top and bottom 64 bit halves.
|
| So base b=2^32 works for C but not Fortran,
| and base B=2^64 works for neither C nor Fortran.
|
| I offer here is a prime that provides CMWC RNGs for both
| 32- and 64-bits, and for both C and Fortran, and with
| equally massive periods, again greater than 2^(1.3million):
|
| p=640*b^41265+1 = 2748779069440*B^20632+1 = 5*2^1320487+1.
|
| That prime came from the many who have dedicated their
| efforts and computer time to prime searches. After some
| three weeks of dedicated computer time using pfgw with
| scrypt, I found the orders of b and B:
| 5*2^1320481 for b=2^32, 5*2^1320480 for B=2^64.
|
| It is the choice of the "a" that makes it feasible to get
| the top and bottom valves of t=a*x+c, yet stay within the
| integer sizes the C or Fortran compilers are set for.
| In the above prime: a=640=2^9+2^7 for b=2^32 and
| a=2748779069440=2^41+2^39 for B=2^64.
| Thus, for example with b=2^32 and using only 32-bit C code,
| with a supposed 128-bit t=(2^9+2^7)*x+c, the top and bottom
| 32-bits of t may be obtained by setting, say,
| h=(c&1); z=(x<<9)>>1 + (x<<7)>>1 + c>>1;
| then the top half of that t would be
| c=(x>>23)+(x>>25)+(z>>31);
| and the bottom half, before being complemented, would be
| x=(z<<1)+h;
|
| When B=2^64 we need only change to
| h=(c&1); z=(x<<41)>>1 + (x<<39)>>1 + c>>1;
| c=(x>>23)+(x>>25)+(z>>63);
|
| These C operations all have Fortran equivalents, and will
| produce the required bit patterns, whether integers are
| considered signed or unsigned. (In C, one must make sure
| that the >> operation performs a logical right shift,
| perhaps best done via "unsigned" declarations.)
|
| The CMWC z "vector" elements [x_1,x_2,...,x_r] are kept in
| an array, Q[] in C, Q() in Fortran, with a separate current
| "carry". This is all spelled out in the following examples:
| code for 32- and 64-bit SuperKiss RNGs for C and Fortran.
|
| Note that in these sample listings, the Q array is seeded
| by CNG+XS, based on the seed values specified in the
| initial declarations. For many simulation studies, the
| 73 bits needed to seed the initial xcng, xs and carry<a
| for the 32-bit version, or 169 bits needed for the 64-bit
| version, may be adequate.
| But more demanding applications may require a significant
| portion of the >1.3 million seed bits that Q requires.
| See text and comments from the Nov 3 posting.
|
| I am indebted to an anonymous mecej4 for providing the basic
| form and KIND declarations of the Fortran versions.
|
| Please let me and other readers know if the results are not
| as specified when run with your compilers, or if you can
| provide equivalent versions in other programming languages.
|
| George Marsaglia
|
| --------------------------------------------------------
| Here is SUPRKISS64.c, the immense-period 64-bit RNG. I
| invite you to cut, paste, compile and run to see if you
| get the result I do. It should take around 20 seconds.
| --------------------------------------------------------
| /* SUPRKISS64.c, period 5*2^1320480*(2^64-1) */
| #include <stdio.h>
| static unsigned long long Q[20632],carry=36243678541LL,
| xcng=12367890123456LL,xs=521288629546311LL,indx=20632;
|
| #define CNG ( xcng=6906969069LL*xcng+123 )
| #define XS ( xs^=xs<<13,xs^=xs>>17,xs^=xs<<43 )
| #define SUPR ( indx<20632 ? Q[indx++] : refill() )
| #define KISS SUPR+CNG+XS
|
| unsigned long long refill( )
| {int i; unsigned long long z,h;
| for(i=0;i<20632;i++){ h=(carry&1);
| z=((Q[i]<<41)>>1)+((Q[i]<<39)>>1)+(carry>>1);
| carry=(Q[i]>>23)+(Q[i]>>25)+(z>>63);
| Q[i]=~((z<<1)+h); }
| indx=1; return (Q[0]);
| }
|
| int main()
| {int i; unsigned long long x;
| for(i=0;i<20632;i++) Q[i]=CNG+XS;
| for(i=0;i<1000000000;i++) x=KISS;
| printf("Does x=4013566000157423768\n x=%LLd.\n",x);
| }
| ---------------------------------------------------------
|
| Here is SUPRKISS32.c, the immense-period 32-bit RNG. I
| invite you to cut, paste, compile and run to see if you
| get the result I do. It should take around 10 seconds.
| ---------------------------------------------------------
| /*suprkiss64.c
| b=2^64; x[n]=(b-1)-[(2^41+2^39)*x[n-20632]+carry mod b]
| period 5*2^1320480>10^397505
| This version of SUPRKISS doesn't use t=a*x+c in 128 bits,
| but uses only 64-bit stuff, takes 20 nanos versus 7.5 for
| the 32-bit unsigned long long t=a*x+c version.
| */
|
| /* SUPRKISS64.c, period 5*2^1320480*(2^64-1) */
| #include <stdio.h>
| static unsigned long long Q[20632],carry=36243678541LL,
| xcng=12367890123456LL,xs=521288629546311LL,indx=20632;
|
| #define CNG ( xcng=6906969069LL*xcng+123 )
| #define XS ( xs^=xs<<13,xs^=xs>>17,xs^=xs<<43 )
| #define SUPR ( indx<20632 ? Q[indx++] : refill() )
| #define KISS SUPR+CNG+XS
|
| unsigned long long refill( )
| {int i; unsigned long long z,h;
| for(i=0;i<20632;i++){ h=(carry&1);
| z=((Q[i]<<41)>>1)+((Q[i]<<39)>>1)+(carry>>1);
| carry=(Q[i]>>23)+(Q[i]>>25)+(z>>63);
| Q[i]=~((z<<1)+h); }
| indx=1; return (Q[0]);
| }
|
| int main()
| {int i; unsigned long long x;
| for(i=0;i<20632;i++) Q[i]=CNG+XS;
| for(i=0;i<1000000000;i++) x=KISS;
| printf("Does x=4013566000157423768\n x=%LLd.\n",x);
| }
|
| -----------------------------------------------------------
|
| And here are equivalent Fortran versions, which, absent
| C's inline features, seem to need ~10% more run time.
|
| -----------------------------------------------------------
| module suprkiss64_M ! period 5*2^1320480*(2^64-1)
| integer,parameter :: I8=selected_int_kind(18)
| integer(kind=I8) :: Q(20632),carry=36243678541_I8, &
| xcng=12367890123456_I8,xs=521288629546311_I8,indx=20633_I8
| contains
| function KISS64() result(x)
| integer(kind=I8) :: x
| if(indx <= 20632)then; x=Q(indx); indx=indx+1
| else; x=refill(); endif
| xcng=xcng*6906969069_I8+123
| xs=ieor(xs,ishft(xs,13))
| xs=ieor(xs,ishft(xs,-17))
| xs=ieor(xs,ishft(xs,43))
| x=x+xcng+xs
| return; end function KISS64
|
| function refill() result(s)
| integer(kind=I8) :: i,s,z,h
| do i=1,20632
| h=iand(carry,1_I8)
| z = ishft(ishft(Q(i),41),-1)+ &
| ishft(ishft(Q(i),39),-1)+ &
| ishft(carry,-1)
| carry=ishft(Q(i),-23)+ishft(Q(i),-25)+ishft(z,-63)
| Q(i)=not(ishft(z,1)+h)
| end do
| indx=2; s=Q(1)
| return; end function refill
|
| end module suprkiss64_M
|
| program testKISS64
| use suprkiss64_M
| integer(kind=I8) :: i,x
| do i=1,20632 !fill Q with Congruential+Xorshift
| xcng=xcng*6906969069_I8+123
| xs=ieor(xs,ishft(xs,13))
| xs=ieor(xs,ishft(xs,-17))
| xs=ieor(xs,ishft(xs,43))
| Q(i)=xcng+xs
| end do
| do i=1,1000000000_I8; x=KISS64(); end do
| write(*,10) x
| 10 format(' Does x = 4013566000157423768 ?',/,6x,'x = ',I20)
| end program testKISS64
| -------------------------------------------------------------
|
| module suprkiss32_M ! period 5*2^1320481*(2^32-1)
| integer,parameter :: I4=selected_int_kind(9)
| integer(kind=I4) :: Q(41265),carry=362_I4, &
| xcng=1236789_I4,xs=521288629_I4,indx=41266_I4
| contains
| function KISS32() result(x)
| integer(kind=I4):: x
| if(indx <= 41265)then;x=Q(indx); indx=indx+1
| else; x=refill(); endif
| xcng=xcng*69069_I4+123
| xs=ieor(xs,ishft(xs,13))
| xs=ieor(xs,ishft(xs,-17));
| xs=ieor(xs,ishft(xs,5))
| x=x+xcng+xs
| return; end function KISS32
|
| function refill() result(s)
| integer(kind=I4) :: i,s,z,h
| do i = 1,41265
| h = iand(carry,1_I4)
| z = ishft(ishft(Q(i),9),-1)+ &
| ishft(ishft(Q(i),7),-1)+ &
| ishft(carry,-1)
| carry=ishft(Q(i),-23)+ishft(Q(i),-25)+ishft(z,-31)
| Q(i)=not(ishft(z,1)+h)
| end do
| indx=2; s=Q(1)
| return; end function refill
|
| end module suprkiss32_M
|
| program testKISS32
| use suprkiss32_M
| integer(kind=I4) :: i,x
| do i=1,41265 !fill Q with Congruential+Xorshift
| xcng=xcng*69069_I4+123
| xs=ieor(xs,ishft(xs,13))
| xs=ieor(xs,ishft(xs,-17))
| xs=ieor(xs,ishft(xs,5))
| Q(i)=xcng+xs
| end do
| do i=1,1000000000_I4; x=KISS32(); end do
| write(*,10) x
| 10 format(' Does x = 1809478889 ?',/,6x,'x =',I11)
| end program testKISS32
|
| ---------------------------------------------------------------
Here is the 32-bit version which I translated:
testKISS32: procedure options (main); /* period 5*2^1320481*(2^32-1) */
declare ( Q(41265), carry initial (362),
xcng initial (1236789), xs initial (521288629),
indx initial (41266)) fixed binary (31);
KISS32: procedure() returns(fixed binary (31)) options (reorder);
declare x fixed binary (31);
if indx <= 41265 then do; x=Q(indx); indx=indx+1; end;
else x=refill();
xcng=xcng*69069+123;
xs=ieor(xs,isll(xs,13));
xs=ieor(xs,isrl(xs,17));
xs=ieor(xs,isll(xs,5));
return (x+xcng+xs);
end KISS32;
refill: procedure() returns(fixed binary(31)) options (reorder);
declare (i,s,z,h) fixed binary (31);
do i = 1 to 41265;
h = iand(carry,1);
z = isrl(isll(Q(i),9),1)+
isrl(isll(Q(i),7),1)+
isrl(carry,1);
carry=isrl(Q(i),23)+isrl(Q(i),25)+isrl(z,31);
Q(i)=inot(isll(z,1)+h);
end;
indx=2;
return (Q(1));
end refill;
declare (i,x) fixed binary (31);
do i=1 to 41265; /* fill Q with Congruential+Xorshift */
xcng=xcng*69069+123;
xs=ieor(xs,isll(xs,13));
xs=ieor(xs,isrl(xs,17));
xs=ieor(xs,isll(xs,5));
Q(i)=xcng+xs;
end;
do i=1 to 1000000000; x=KISS32(); end;
put skip edit ( 'Does x = 1809478889 ?', 'x =', x)
(a, skip, x(5), a, f(11));
end testKISS32;
/*--------------------------------------------------------------- */
robin
| On Nov 3 I posted
|
| RNGs: A Super KISS sci.math, comp.lang.c, sci.crypt
|
| a KISS (Keep-It-Simple-Stupid) RNG combining,
| by addition mod 2^32, three simple RNGs:
| CMWC(Complementary-Multiply-With-Carry)+CNG(Congruential)
| +XS(Xorshift)
| with resulting period greater than 10^402575.
Here is the 64-bit PL/I version which I translated:
testKISS64: procedure options (main); /* period 5*2^1320480*(2^64-1) */
declare (i,x) fixed binary (63);
declare (Q(20632), carry initial (36243678541),
xcng initial (12367890123456), xs initial (521288629546311),
indx initial (20633)) fixed binary (63);
KISS64: procedure() returns (fixed binary (63)) options (reorder);
declare x fixed binary (63);
if indx <= 20632 then do; x=Q(indx); indx=indx+1; end;
else x=refill();
xcng=xcng*6906969069+123;
xs=ieor(xs,isll(xs,13));
xs=ieor(xs,isrl(xs,17));
xs=ieor(xs,isll(xs,43));
return (x+xcng+xs);
end KISS64;
refill: procedure () returns (fixed binary (63)) options (reorder);
declare (i,s,z,h) fixed binary (63);
do i=1 to 20632;
h=iand(carry,1);
z = isrl(isll(Q(i),41),1)+
isrl(isll(Q(i),39),1)+
isrl(carry,1);
carry=isrl(Q(i),23)+isrl(Q(i),25)+isrl(z,63);
Q(i)=inot(isll(z,1)+h);
end;
indx=2;
return(Q(1));
end refill;
do i=1 to 20632; /* fill Q with Congruential+Xorshift */
xcng=xcng*6906969069+123;
xs=ieor(xs,isll(xs,13));
xs=ieor(xs,isrl(xs,17));
xs=ieor(xs,isll(xs,43));
Q(i)=xcng+xs;
end;
do i=1 to 1000000000; x=KISS64(); end;
put skip edit ( 'Does x = 4013566000157423768 ?', 'x = ', x)
(a, skip, x(5), a, f(20));
end testKISS64;
/*-------------------------------------------------------------*/
robin
/*-------------------------------------------------------------*/
testKISS32: procedure options (main);
/* George Marsaglia's Pseudo-Random Number Generator having a period of */
/* 5*2^1320481*(2^32-1). */
/* 32-bit PRNGs are generated directly. */
/* The procedure FILL must be called first before invoking KISS32. */
KISS32: procedure() returns(fixed binary (32) unsigned) options (reorder);
declare x fixed binary (32) unsigned;
declare Q(41265) fixed binary (32) unsigned external,
indx fixed binary (31) static initial (41266);
declare (carry initial (362),
xcng initial (1236789),
xs initial (521288629)) fixed binary (32) unsigned external;
if indx <= 41265 then do; x=Q(indx); indx=indx+1; end;
else do; x=refill(); indx = 2; end;
xcng=xcng*69069+123;
xs=ieor(xs,isll(xs,13));
xs=ieor(xs,isrl(xs,17));
xs=ieor(xs,isll(xs,5));
return (x+xcng+xs);
end KISS32;
refill: procedure() returns(fixed binary(32) unsigned) options (reorder);
declare (i,s,z,h) fixed binary (32) unsigned;
declare Q(41265) fixed binary (32) unsigned external;
declare (carry initial (362),
xcng initial (1236789),
xs initial (521288629)) fixed binary (32) unsigned external;
do i = 1 to 41265;
h = iand(carry,1);
z = isrl(isll(Q(i),9),1) + isrl(isll(Q(i),7),1) + isrl(carry,1);
carry=isrl(Q(i),23)+isrl(Q(i),25)+isrl(z,31);
Q(i)=inot(isll(z,1)+h);
end;
return (Q(1));
end refill;
fill: procedure options (reorder);
declare Q(41265) fixed binary (32) unsigned external;
declare (carry initial (362),
xcng initial (1236789),
xs initial (521288629)) fixed binary (32) unsigned external;
declare i fixed binary (31);
do i=1 to 41265; /* fill Q with Congruential+Xorshift */
xcng=xcng*69069+123;
xs=ieor(xs,isll(xs,13));
xs=ieor(xs,isrl(xs,17));
xs=ieor(xs,isll(xs,5));
Q(i)=xcng+xs;
end;
end fill;
declare i fixed binary (31), x fixed binary (32) unsigned;
call fill;
robin
testKISS64: procedure options (main);
declare (i,x) fixed binary (64) unsigned;
/* George Marsaglia's Pseudo-Random Number Generator having a period of */
/* 5*2^1320480*(2^64-1). */
/* 64-bit PRNGs are generated directly. */
/* The procedure FILL must be called first before invoking KISS64. */
KISS64: procedure() returns (fixed binary (64) unsigned) options (reorder);
declare x fixed binary (64) unsigned,
Q(20632) fixed binary (64) unsigned external,
indx fixed binary static initial (20633);
declare (carry initial (36243678541),
xcng initial (12367890123456),
xs initial (521288629546311)) fixed binary (64) unsigned external;
if indx <= 20632 then do; x=Q(indx); indx=indx+1; end;
else do; x=refill(); indx = 2; end;
(nofixedoverflow):
xcng=xcng*6906969069+123;
xs=ieor(xs,isll(xs,13));
xs=ieor(xs,isrl(xs,17));
xs=ieor(xs,isll(xs,43));
return (x+xcng+xs);
end KISS64;
refill: procedure () returns (fixed binary (64) unsigned) options (reorder);
declare Q(20632) fixed binary (64) unsigned external;
declare (carry initial (36243678541),
xcng initial (12367890123456),
xs initial (521288629546311)) fixed binary (64) unsigned external;
declare i fixed binary, (z,h) fixed binary (64) unsigned;
do i=1 to 20632;
h=iand(carry,1);
z = isrl(isll(Q(i),41),1) + isrl(isll(Q(i),39),1) + isrl(carry,1);
carry=isrl(Q(i),23)+isrl(Q(i),25)+isrl(z,63);
Q(i)=inot(isll(z,1)+h);
end;
return (Q(1));
end refill;
fill: procedure options (reorder);
declare Q(20632) fixed binary (64) unsigned external;
declare (carry initial (36243678541),
xcng initial (12367890123456),
xs initial (521288629546311)) fixed binary (64) unsigned external;
declare i fixed binary;
do i=1 to 20632; /* fill Q with Congruential+Xorshift */
xcng=xcng*6906969069+123;
xs=ieor(xs,isll(xs,13));
xs=ieor(xs,isrl(xs,17));
xs=ieor(xs,isll(xs,43));
Q(i)=xcng+xs;
end;
end fill;
call fill;
do i=1 to 1000000000; x=KISS64(); end;
put skip edit ( 'Does x = 4013566000157423768 ?', 'x = ', x)
(a, skip, x(5), a, f(20));
end testKISS64;
/*-------------------------------------------------------------*/
robin
/*-----------------------------------------------------------*/
RNG64: package exports (KISS64, FILL);
declare Q(20632) fixed binary (64) unsigned static;
declare (carry initial (36243678541),
xcng initial (12367890123456),
xs initial (521288629546311)) fixed binary (64) unsigned static;
/* George Marsaglia's Pseudo-Random Number Generator having a period of */
/* 5*2^1320480*(2^64-1). */
/* 64-bit PRNGs are generated directly. */
/* The procedure FILL must be called first before invoking KISS64. */
/* Translated from Fortran to PL/I by R. Vowels. 29 November 2009. */
KISS64: procedure() returns (fixed binary (64) unsigned) options (reorder);
declare x fixed binary (64) unsigned,
indx fixed binary static initial (20633);
if indx <= 20632 then do; x=Q(indx); indx=indx+1; end;
else do; x=refill(); indx = 2; end;
(nofixedoverflow):
xcng=xcng*6906969069+123;
xs=ieor(xs,isll(xs,13));
xs=ieor(xs,isrl(xs,17));
xs=ieor(xs,isll(xs,43));
return (x+xcng+xs);
end KISS64;
refill: procedure () returns (fixed binary (64) unsigned) options (reorder);
declare i fixed binary, (z,h) fixed binary (64) unsigned;
do i=1 to 20632;
h=iand(carry,1);
z = isrl(isll(Q(i),41),1) + isrl(isll(Q(i),39),1) + isrl(carry,1);
carry=isrl(Q(i),23)+isrl(Q(i),25)+isrl(z,63);
Q(i)=inot(isll(z,1)+h);
end;
return (Q(1));
end refill;
fill: procedure options (reorder);
declare i fixed binary;
do i=1 to 20632; /* fill Q with Congruential+Xorshift */
xcng=xcng*6906969069+123;
xs=ieor(xs,isll(xs,13));
xs=ieor(xs,isrl(xs,17));
xs=ieor(xs,isll(xs,43));
Q(i)=xcng+xs;
end;
end fill;
end RNG64;
/*-------------------------------------------------------------*/