1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
|
/* This program by D E Knuth is in the public domain and freely copyable.
* It is explained in Seminumerical Algorithms, 3rd edition, Section 3.6
* (or in the errata to the 2nd edition --- see
* http://www-cs-faculty.stanford.edu/~knuth/taocp.html
* in the changes to Volume 2 on pages 171 and following). */
/* N.B. The MODIFICATIONS introduced in the 9th printing (2002) are
included here; there's no backwards compatibility with the original. */
/* This version also adopts Brendan McKay's suggestion to
accommodate naive users who forget to call ran_start(seed). */
/* If you find any bugs, please report them immediately to
* [email protected]
* (and you will be rewarded if the bug is genuine). Thanks! */
/************ see the book for explanations and caveats! *******************/
/************ in particular, you need two's complement arithmetic **********/
#define KK 100 /* the long lag */
#define LL 37 /* the short lag */
#define MM (1L<<30) /* the modulus */
#define mod_diff(x,y) (((x)-(y))&(MM-1)) /* subtraction mod MM */
long ran_x[KK]; /* the generator state */
#ifdef __STDC__
void ran_array(long aa[],int n)
#else
void ran_array(aa,n) /* put n new random numbers in aa */
long *aa; /* destination */
int n; /* array length (must be at least KK) */
#endif
{
register int i,j;
for (j=0;j<KK;j++) aa[j]=ran_x[j];
for (;j<n;j++) aa[j]=mod_diff(aa[j-KK],aa[j-LL]);
for (i=0;i<LL;i++,j++) ran_x[i]=mod_diff(aa[j-KK],aa[j-LL]);
for (;i<KK;i++,j++) ran_x[i]=mod_diff(aa[j-KK],ran_x[i-LL]);
}
/* the following routines are from exercise 3.6--15 */
/* after calling ran_start, get new randoms by, e.g., "x=ran_arr_next()" */
#define QUALITY 1009 /* recommended quality level for high-res use */
long ran_arr_buf[QUALITY];
long ran_arr_dummy=-1, ran_arr_started=-1;
long *ran_arr_ptr=&ran_arr_dummy; /* the next random number, or -1 */
#define TT 70 /* guaranteed separation between streams */
#define is_odd(x) ((x)&1) /* units bit of x */
#ifdef __STDC__
void ran_start(long seed)
#else
void ran_start(seed) /* do this before using ran_array */
long seed; /* selector for different streams */
#endif
{
register int t,j;
long x[KK+KK-1]; /* the preparation buffer */
register long ss=(seed+2)&(MM-2);
for (j=0;j<KK;j++) {
x[j]=ss; /* bootstrap the buffer */
ss<<=1; if (ss>=MM) ss-=MM-2; /* cyclic shift 29 bits */
}
x[1]++; /* make x[1] (and only x[1]) odd */
for (ss=seed&(MM-1),t=TT-1; t; ) {
for (j=KK-1;j>0;j--) x[j+j]=x[j], x[j+j-1]=0; /* "square" */
for (j=KK+KK-2;j>=KK;j--)
x[j-(KK-LL)]=mod_diff(x[j-(KK-LL)],x[j]),
x[j-KK]=mod_diff(x[j-KK],x[j]);
if (is_odd(ss)) { /* "multiply by z" */
for (j=KK;j>0;j--) x[j]=x[j-1];
x[0]=x[KK]; /* shift the buffer cyclically */
x[LL]=mod_diff(x[LL],x[KK]);
}
if (ss) ss>>=1; else t--;
}
for (j=0;j<LL;j++) ran_x[j+KK-LL]=x[j];
for (;j<KK;j++) ran_x[j-LL]=x[j];
for (j=0;j<10;j++) ran_array(x,KK+KK-1); /* warm things up */
ran_arr_ptr=&ran_arr_started;
}
#define ran_arr_next() (*ran_arr_ptr>=0? *ran_arr_ptr++: ran_arr_cycle())
long ran_arr_cycle()
{
if (ran_arr_ptr==&ran_arr_dummy)
ran_start(314159L); /* the user forgot to initialize */
ran_array(ran_arr_buf,QUALITY);
ran_arr_buf[KK]=-1;
ran_arr_ptr=ran_arr_buf+1;
return ran_arr_buf[0];
}
#include <stdio.h>
int main()
{
register int m; long a[2009];
ran_start(310952L);
for (m=0;m<=2009;m++) ran_array(a,1009);
printf("%ld\n", a[0]); /* 995235265 */
ran_start(310952L);
for (m=0;m<=1009;m++) ran_array(a,2009);
printf("%ld\n", a[0]); /* 995235265 */
return 0;
}
|
