Rev 237 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
| Rev | Author | Line No. | Line |
|---|---|---|---|
| 259 | daniel-mar | 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 | } |