How to seed and use SUPRKISS64.c?

[kfmfe04]

by the venerable Professor George Marsaglia who passed away just 4 months ago

source:
http://forums.silverfrost.com/viewtopic.php?t=1480

question: how do you seed and use this thing like srand_sk64() and rand_sk64()?
kev82? anyone else?

changes:
- formatting
- line 30 LL format changed to ll, output changed to 4013566000157423768, as per thread

with the fix:
- line 15 changed to 20632 (from 20635) to prevent overflow

Does 4013566000157423768
 x=4013566000157423768.

/* 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); 
}

[kev82]

As far as I can tell, it is made up of three generators.

The first is seeded by the variable xchg, and I've checked it has full period for all seeds, so you can seed it with any 64bit value.

I'm not familiar with the second, but it is seeded by xs. I assume you can seed it with any value - because there should be an attached warning if you couldn't.

The third is seeded by the for loop on line 28, so you needn't worry about that.

[kfmfe04]

kev82, thx for your quick reply - I see that I should leave Q and indx alone

in this thread:

http://www.velocityreviews.com/forums/t706893-superkiss-for-32-and-64-bit-rngs-in-both-c-and-fortran.html

it mentions that the 64-bit version uses 169-bits (64 + 73 + 32?)

so out of the 3 variables left, carry, xcng, and xs, you think I should leave carry and xs alone and just seed xcng?

in the thread above, note:

gmarsaglia wrote:

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.

edit: also, I assume that every time I call KISS, it's the equivalent of calling rand_sk64(), correct?

[kev82]

I can't figure the 73 (in the 32bit version) or 169 in this one. It may be worth posting on that forum or one of the newsgroups for an explanation. As I said, I don't work on this area any more.

I wouldn't alter index as that is set to force a call to refill on the first invocation of XS.

It is the `a' that is confusing me. I understand a to be 6906969069, but when you then take 64+64+log2(6906969069) you don't get 169.

Yes, KISS is the macro that does the generation.

[kfmfe04]

I am trying to understand this posting now:

http://www.adras.com/RNGs-A-Super-KISS.t567171-93-2.html

btw, if you don't mind me asking, is your primary background in math, CS, or both? what specialties?

edit: I'm beginning to think I need to seed xcng and xs, where xs<a (this should set up Q, I think)

[kev82]

Good luck, I'm afraid it's a little beyond me.

My background is maths. I have no formal CS training of any kind. My masters thesis was on distribution transforms. That's basically once you've made the effort to get a uniform number from a prng what do you then do if you want a normal/exponential/weibull/etc number.

For the last few years I have been working on a particular type of VRP using metaheuristic approaches.

[kfmfe04]

that explains your aptitude in this area!

ok - that last posting is getting me a little bit closer into figuring out how to set this up

the actual seed appears to be all the bits in Q[] (20632 * 64 = 1.3mln) and the initial carry, but if one doesn't need the extra power, one could bootstrap Q[] with xcng and xs as shown on Line 28 - I think the 64 bits George refers to below in xcng and xs is 2 x 32 bit in the 32 bit version.

edit: 1337280 / 64 = 20895 (wonder where this 20895 came from)

gmarsaglia wrote:

"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.

[kev82]

the actual seed appears to be all the bits in Q[]

Bear in mind the actual generator is CNG + XS + SUPR. The Q array is only used for the calculation of SUPR. So even if you didn't initialize it from CNG and XS, their initial values do play a partin the overall generator.