Rev 237 | Go to most recent revision | Show entire file | Regard whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 237 | Rev 259 | ||
---|---|---|---|
Line 1... | Line -... | ||
1 | /* This program by D E Knuth is in the public domain and freely copyable |
- | |
2 | * AS LONG AS YOU MAKE ABSOLUTELY NO CHANGES! |
1 | /* This program by D E Knuth is in the public domain and freely copyable |
- | 2 | * AS LONG AS YOU MAKE ABSOLUTELY NO CHANGES! |
|
- | 3 | * It is explained in Seminumerical Algorithms, 3rd edition, Section 3.6 |
|
- | 4 | * (or in the errata to the 2nd edition, see |
|
- | 5 | * https://www-cs-faculty.stanford.edu/~knuth/taocp.html |
|
- | 6 | * in the changes to pages 171 and following). |
|
- | 7 | * |
|
- | 8 | * If you find any bugs, please report them immediately to |
|
- | 9 | * taocp@cs.stanford.edu |
|
- | 10 | * (and you will be rewarded if the bug is genuine). Thanks! */ |
|
- | 11 | ||
- | 12 | #define KK 100 /* the long lag */ |
|
- | 13 | #define LL 37 /* the short lag */ |
|
- | 14 | #define MM (1<<30) /* the modulus */ |
|
- | 15 | #define mod_diff(x,y) (x-y)&(MM-1) /* subtraction mod MM */ |
|
- | 16 | ||
- | 17 | long ran_x[KK]; /* the generator state */ |
|
- | 18 | ||
- | 19 | void ran_array(aa,n) /* put n new random numbers in aa */ |
|
- | 20 | long *aa; /* destination */ |
|
- | 21 | int n; /* array length (must be at least KK) */ |
|
- | 22 | { |
|
- | 23 | register int i,j; |
|
- | 24 | for (j=0;j<KK;j++) aa[j]=ran_x[j]; |
|
- | 25 | for (;j<n;j++) aa[j]=mod_diff(aa[j-KK],aa[j-LL]); |
|
- | 26 | for (i=0;i<LL;i++,j++) ran_x[i]=mod_diff(aa[j-KK],aa[j-LL]); |
|
- | 27 | for (;i<KK;i++,j++) ran_x[i]=mod_diff(aa[j-KK],ran_x[i-LL]); |
|
- | 28 | } |
|
- | 29 | ||
- | 30 | #define TT 70 /* guaranteed separation between streams */ |
|
- | 31 | #define is_odd(x) (x&1) /* units bit of x */ |
|
- | 32 | #define evenize(x) ((x)&(MM-2)) /* make x even */ |
|
- | 33 | ||
- | 34 | void ran_start(seed) /* do this before using ran_array */ |
|
- | 35 | long seed; /* selector for different streams */ |
|
- | 36 | { |
|
- | 37 | register int t,j; |
|
- | 38 | long x[KK+KK-1]; /* the preparation buffer */ |
|
- | 39 | register long ss=evenize(seed+2); |
|
- | 40 | for (j=0;j<KK;j++) { |
|
- | 41 | x[j]=ss; /* bootstrap the buffer */ |
|
- | 42 | ss<<=1; if (ss>=MM) ss-=MM-2; /* cyclic shift 29 bits */ |
|
- | 43 | } |
|
- | 44 | for (;j<KK+KK-1;j++) x[j]=0; |
|
- | 45 | x[1]++; /* make x[1] (and only x[1]) odd */ |
|
- | 46 | ss=seed; |
|
- | 47 | t=TT-1; while (t) { |
|
- | 48 | for (j=KK-1;j>0;j--) x[j+j]=x[j]; /* "square" */ |
|
- | 49 | for (j=KK+KK-2;j>KK-LL;j-=2) x[KK+KK-1-j]=evenize(x[j]); |
|
- | 50 | for (j=KK+KK-2;j>=KK;j--) if(is_odd(x[j])) { |
|
- | 51 | x[j-(KK-LL)]=mod_diff(x[j-(KK-LL)],x[j]); |
|
- | 52 | x[j-KK]=mod_diff(x[j-KK],x[j]); |
|
- | 53 | } |
|
- | 54 | if (is_odd(ss)) { /* "multiply by z" */ |
|
- | 55 | for (j=KK;j>0;j--) x[j]=x[j-1]; |
|
- | 56 | x[0]=x[KK]; /* shift the buffer cyclically */ |
|
- | 57 | if (is_odd(x[KK])) x[LL]=mod_diff(x[LL],x[KK]); |
|
- | 58 | } |
|
- | 59 | if (ss) ss>>=1; else t--; |
|
- | 60 | } |
|
- | 61 | for (j=0;j<LL;j++) ran_x[j+KK-LL]=x[j]; |
|
- | 62 | for (;j<KK;j++) ran_x[j-LL]=x[j]; |
|
- | 63 | } |
|
- | 64 |