|
|
1.1 root 1: #define P 98
2: #define Q 27
3: #define DELAY 9800 /* 98*40 might be adequate */
4:
5: /* Feedback shift register random number generator.
6: * Steve Zucker -- March 5, 1975
7: * Period is 2**98 - 1.
8: * See JACM, Vol 20 No 3, July, 1973, pp. 456-468.
9: * Initialization modified: J. Gillogly, 11 Mar 75
10: * Modified to use long ints 12 Apr 77
11: */
12:
13: int rand_t[P],rand_x[P];
14: int *rand_end = {&rand_t[P]};
15: int *rand_j = {&rand_t[P]};
16:
17: random()
18: { register int *j, *k, *l;
19: if ( (j = ++rand_j) >= (k = rand_end) )
20: rand_j = j = rand_t;
21: return (*j ^= ( (l = j+Q) >= k ? *(l-P) : *l) );
22: }
23:
24:
25: #define TABSIZE 8000 /* table of initial values */
26: char *ivp; /* pointer to the table */
27: int pshift; /* will be set to P/16 */
28: int mask[16]; /* mask[0] has a 1-bit in the 0 position, etc */
29: long p2n[27], pq2n[27];
30:
31: srandom(n) /* compatible with unix initialization */
32: int n;
33: { char buf[14];
34: sprintf(buf,"%d",n & 0xffff); /* convert number to string */
35: initr(buf); /* and use that to initialize */
36: }
37:
38: initr(c)
39: char *c;
40: { register int col, *k, j;
41: long offset;
42: int m,l,*i,deb;
43: auto char iv[TABSIZE];
44: auto int column[P/16+1];
45: pshift=P/16;
46: for (j=0; j<16; j++) mask[j]=1<<j;
47: p2n[0]=P; /* store 2**N*P */
48: pq2n[0]=P-Q; /* and 2**N*(P-Q) */
49: for (j=1; j<27; j++)
50: { p2n[j]=p2n[j-1]<<1;
51: pq2n[j]=pq2n[j-1]<<1;
52: }
53: if (c)
54: { for (i=rand_x, j=0, ivp=iv; j<P; j++) /* bits in rand_x */
55: *ivp++ = *i++ = (*(c+(j>>3)) >> (j&07)) & 01;
56: for (k=rand_t; k<rand_end;) *k++ = 0;
57: for (ivp = &iv[P]; ivp<&iv[TABSIZE]; ivp++) /* init table */
58: *ivp = *(ivp-P+Q) ^ *(ivp-P);
59: ivp=iv; /* point to table so find can find it */
60: offset=P*5; /* getting up to P*5000 */
61: offset=(offset<<10)-P*120; /* NG for higher offsets */
62: for (l=15; l--;) /* for each column */
63: { find(offset+=DELAY,column); /* result stored in column */
64: for (j=0; j<P; j++)
65: { col=column[j>>4]; /* pick right result*/
66: m = 1<<(15-(j%16)); /* and bit posn */
67: if ((col&m) != 0) /* the bit is on */
68: rand_t[j] |= 1<<l; /* turn on in init */
69: }
70: }
71: for (k=rand_x, i=rand_t; i<rand_end;) *k++ = *i++;
72: }
73: rand_j = rand_end;
74: for (i=rand_t, k=rand_x; i<rand_end;) *i++ = *k++;
75: }
76:
77: find(offset,result)
78: long offset;
79: int *result; /* P/16-long column of results */
80: { register int n,l;
81: register char *i;
82: long j,k;
83: auto int lcol[P/16+1],rcol[P/16+1]; /* to save results in */
84: if (offset<TABSIZE-P)
85: { i=ivp+(n=offset); /* pt at right place in iv */
86: for (l=0; l<=pshift; l++) /* abt P/16 result words */
87: { result[l]=0; /* turn off all bits */
88: for (n=0; n<16; n++) /* cycle thru the bits in wd*/
89: if (*i++) /* next initial bit on */
90: result[l] |= mask[15-n];
91: /* put next initial bit into the right result bit */
92: }
93: return;
94: }
95: /* now find max n for which offset-p*2**n>0, and will use that n */
96: for (n=0; p2n[n]<=offset; n++)
97: if (n>26) printf("ERROR: exceeded 2**N*P table\n");
98: if (n==0) printf("ERROR: finding max after done!\n");
99: j=offset-pq2n[n-1];
100: k=offset-p2n[n-1];
101: if (j<TABSIZE-P) /* no need to call again */
102: { i=ivp+(n=j); /* pt at right place in iv */
103: for (l=0; l<=pshift; l++) /* abt P/16 result words */
104: { lcol[l]=0; /* turn off all bits */
105: for (n=0; n<16; n++) /* cycle thru the bits in wd*/
106: if (*i++) /* next initial bit on */
107: lcol[l] |= mask[15-n];
108: /* put next initial bit into the right result bit */
109: }
110: }
111: else find(j,lcol);
112: if (k<TABSIZE-P) /* no need to call again */
113: { i=ivp+(n=k); /* pt at right place in iv */
114: for (l=0; l<=pshift; l++) /* abt P/16 result words */
115: { rcol[l]=0; /* turn off all bits */
116: for (n=0; n<16; n++) /* cycle thru the bits in wd*/
117: if (*i++) /* next initial bit on */
118: rcol[l] |= mask[15-n];
119: /* put next initial bit into the right result bit */
120: }
121: }
122: else find(k,rcol);
123: for (l=0; l<=P>>4; l++) result[l] = lcol[l]^rcol[l];
124: }
125:
126:
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.