Subversion Repositories filter_foundry

Rev

Rev 237 | Go to most recent revision | Blame | Last modification | View Log | RSS feed

  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.