Prev: Positive definite and positive definite submatrices
Next: Why G. Chaitin mesure of complexity is wrong.
From: geo on 27 Nov 2009 16:52
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

---------------------------------------------------------------
From: Ben Bacarisse on 27 Nov 2009 17:17
geo <gmarsaglia(a)gmail.com> writes:
<snip interesting stuff>
> --------------------------------------------------------
> 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.
> --------------------------------------------------------
<snip>
> 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);

Standard C uses %lld for printing long long int. Changing %LLd to
%llu (u is correct for unsigned types like x) I get the expected
output from gcc 4.4.1 on 64-bit Intel Linux.

> }

<snip more>
--
Ben.
From: robin on 28 Nov 2009 17:10
"geo" <gmarsaglia(a)gmail.com> wrote in message news:d32ec0cc-5937-4082-8322-0111711307c4(a)j24g2000yqa.googlegroups.com...
| 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.

Fantastic!



From: amzoti on 29 Nov 2009 22:35
On Nov 27, 1:52 pm, geo <gmarsag...(a)gmail.com> wrote:
>  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))
> ...
>
> read more »

I am a big fan of your work.

One question, has the crypto usage of this been statistically verified
with http://www.iro.umontreal.ca/~simardr/testu01/tu01.html?

I supposed you already tested it against KISS.

Thanks ~A
From: Richard Outerbridge on 29 Nov 2009 23:56
Here are two somewhat "whiter" and more C99 compliant versions of
SuperKISS32 and SuperKISS64.

On my Macintosh PPC PowerBook5,8 running Leopard v10.5.8 (9L30) PPC,
and using for a compiler powerpc-apple-darwin9-gcc-4.0.1, I get UNIX
time speeds of ~11.0 and ~34.9 seconds, respectively. I've assumed
the standard C type restrictions: an int is only at least 16-bits,
and long is only at least 32-bits; and so on.

/* SUPRKISS32.c, period 5*2^1320481*(2^32-1) */

#include <stdlib.h>
#include <stdio.h>

static unsigned long Q[41265], carry = 362UL, \
xcng = 1236789UL, xs = 521288629UL;

static long indx = 41265L;

#define CNG ( xcng = 69069UL * xcng + 123UL )
#define XS ( xs ^= xs << 13, xs ^= xs >> 17, xs ^= xs << 5 )
#define SUPR ( (indx < 41265L) ? Q[indx++] : refill() )
#define KISS SUPR+CNG+XS

unsigned long refill( void ) {
long i;
unsigned long z, h;

for( i = 0; i < 41265L; i++ ) {
h = carry &1UL;
z = ((Q[i] << 9) >> 1) + ((Q[i] << 7) >> 1) + (carry >> 1);
carry = (Q[i] >> 23) + (Q[i] >> 25) + (z >> 31);
Q[i] = ~((z << 1) + h);
}
indx = 1L;
return ( Q[0] );
}

int main( void ) {
long i;
unsigned long x;

for( i = 0; i < 41265L; i++ ) Q[i] = CNG + XS;
for( i = 0; i < 1000000000L; i++ ) x = KISS;
printf("Does x = 1809478889?\n x = %lu.\n", x);

exit(0);
}

/* and */

/* SUPRKISS64.c, period 5*2^1320480*(2^64-1) */

#include <stdlib.h>
#include <stdio.h>

static unsigned long long Q[20632], carry = 36243678541ULL, \
xcng = 12367890123456ULL, xs = 521288629546311ULL;

static int indx = 20632;

#define CNG ( xcng = 6906969069ULL * xcng + 123ULL )
#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( void ) {
int i;
unsigned long long z, h;

for( i = 0; i < 20632; i++ ) {
h = carry &1ULL;
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( void ) {
int i;
long j;
unsigned long long x;

for( i = 0; i < 20632; i++ ) Q[i] = CNG + XS;
for( j = 0; j < 1000000000L; j++ ) x = KISS;
printf("Does x = 4013566000157423768?\n x = %llu.\n", x);

exit(0);
}

Richard Outerbridge
outer(a)interlog.com

--

--�
"Picking random numbers is far too�important�to be left to chance."
 |  Next  |  Last
Pages: 1 2
Prev: Positive definite and positive definite submatrices
Next: Why G. Chaitin mesure of complexity is wrong.