V10/games/atc/ran.c

#define P       98
#define Q       27
#define DELAY   9800    /* 98*40 might be adequate */

/* Feedback shift register random number generator.
 * Steve Zucker -- March 5, 1975
 * Period is 2**98 - 1.
 * See JACM, Vol 20 No 3, July, 1973, pp. 456-468.
 * Initialization modified: J. Gillogly, 11 Mar 75
 * Modified to use long ints 12 Apr 77
 */

int rand_t[P],rand_x[P];
int *rand_end = {&rand_t[P]};
int *rand_j   = {&rand_t[P]};

random()
{       register int *j, *k, *l;
        if ( (j = ++rand_j) >= (k = rand_end) )
                rand_j = j = rand_t;
        return (*j ^= ( (l = j+Q) >= k ? *(l-P) : *l) );
}


#define TABSIZE  8000    /* table of initial values */
char *ivp;              /* pointer to the table */
int pshift;             /* will be set to P/16 */
int mask[16];           /* mask[0] has a 1-bit in the 0 position, etc */
long p2n[27], pq2n[27];

srandom(n)        /* compatible with unix initialization */
int n;
{       char buf[14];
	sprintf(buf,"%d",n & 0xffff);    /* convert number to string */
	initr(buf);             /* and use that to initialize */
}

initr(c)
char *c;
{       register int col, *k, j;
	long offset;
        int m,l,*i,deb;
        auto char iv[TABSIZE];
        auto int column[P/16+1];
        pshift=P/16;
        for (j=0; j<16; j++) mask[j]=1<<j;
	p2n[0]=P;      /* store 2**N*P */
	pq2n[0]=P-Q;  /* and 2**N*(P-Q) */
	for (j=1; j<27; j++)
	{       p2n[j]=p2n[j-1]<<1;
		pq2n[j]=pq2n[j-1]<<1;
	}
        if (c)
        {       for (i=rand_x, j=0, ivp=iv; j<P; j++)   /* bits in rand_x */
                        *ivp++ = *i++ = (*(c+(j>>3)) >> (j&07)) & 01;
                for (k=rand_t; k<rand_end;) *k++ = 0;
                for (ivp = &iv[P]; ivp<&iv[TABSIZE]; ivp++) /* init table */
                        *ivp = *(ivp-P+Q) ^ *(ivp-P);
                ivp=iv; /* point to table so find can find it */
		offset=P*5;     /* getting up to P*5000 */
		offset=(offset<<10)-P*120; /* NG for higher offsets */
                for (l=15; l--;)        /* for each column */
		{       find(offset+=DELAY,column);    /* result stored in column */
                        for (j=0; j<P; j++)
                        {       col=column[j>>4];       /* pick right result*/
                                m = 1<<(15-(j%16));     /* and bit posn */
                                if ((col&m) != 0)       /* the bit is on */
                                    rand_t[j] |= 1<<l;  /* turn on in init */
                        }
                }
                for (k=rand_x, i=rand_t; i<rand_end;) *k++ = *i++;
        }
        rand_j = rand_end;
        for (i=rand_t, k=rand_x; i<rand_end;) *i++ = *k++;
}

find(offset,result)
long offset;
int *result;           /* P/16-long column of results */
{       register int n,l;
        register char *i;
	long j,k;
        auto int lcol[P/16+1],rcol[P/16+1];      /* to save results in */
	if (offset<TABSIZE-P)
        {       i=ivp+(n=offset);               /* pt at right place in iv */
                for (l=0; l<=pshift; l++)       /* abt P/16 result words */
                {       result[l]=0;            /* turn off all bits */
                        for (n=0; n<16; n++)    /* cycle thru the bits in wd*/
                                if (*i++)       /* next initial bit on */
                                        result[l] |= mask[15-n];
                        /* put next initial bit into the right result bit */
                }
                return;
        }
        /* now find max n for which offset-p*2**n>0, and will use that n */
	for (n=0; p2n[n]<=offset; n++)
		if (n>26) printf("ERROR: exceeded 2**N*P table\n");
	if (n==0) printf("ERROR: finding max after done!\n");
	j=offset-pq2n[n-1];
	k=offset-p2n[n-1];
	if (j<TABSIZE-P)                        /* no need to call again */
        {       i=ivp+(n=j);                    /* pt at right place in iv */
                for (l=0; l<=pshift; l++)       /* abt P/16 result words */
                {       lcol[l]=0;              /* turn off all bits */
                        for (n=0; n<16; n++)    /* cycle thru the bits in wd*/
                                if (*i++)       /* next initial bit on */
                                        lcol[l] |= mask[15-n];
                        /* put next initial bit into the right result bit */
                }
        }
        else find(j,lcol);
	if (k<TABSIZE-P)                        /* no need to call again */
        {       i=ivp+(n=k);                    /* pt at right place in iv */
                for (l=0; l<=pshift; l++)       /* abt P/16 result words */
                {       rcol[l]=0;              /* turn off all bits */
                        for (n=0; n<16; n++)    /* cycle thru the bits in wd*/
                                if (*i++)       /* next initial bit on */
                                        rcol[l] |= mask[15-n];
                        /* put next initial bit into the right result bit */
                }
        }
        else find(k,rcol);
        for (l=0; l<=P>>4; l++) result[l] = lcol[l]^rcol[l];
}