From schutz@wehi.edu.au  Sun Mar  3 03:17:10 2002
From: schutz@wehi.edu.au (schutz@wehi.edu.au)
Date: Sun, 3 Mar 2002 04:17:10 +0100 (MET)
Subject: [Rd] Weakness in Knuth-TAOCP RNG (fwd) (PR#1336)
Message-ID: <200203030317.EAA20430@pubhealth.ku.dk>

Prof Brian D Ripley wrote:

>Attachments don't really work with R-bugs.  To save us trying to unpick
>this manually, could you re-send it as the body of a message to R-bugs
>in reply to this one, or quoting (PR#1336) in the subject, please?

Sorry about that (Pine sometimes modify lines of the body without warning,
so I preferred using an attachment). Here is the patch for src/main/RNG.c.

Frdric
---
>> Don Knuth recently updated its RNG, following the discovery of a
>> weakness in the seed initialisation (see
>> http://Sunburn.Stanford.EDU/~knuth/news02.html), so I believe that the
>> corresponding code in R should be updated as well. If it may be usefull,
>> I attach a patch to do the change.
>>
>> Notes:
>>
>>   - only a few lines of code have really changed, but the order of the
>>     functions has changed in Knuth's original code, so the patch is
>>     longer than it should really be;
>>   - Despite Knuth's request that no change should be made to the code,
>>     the present R code (1.5.0-devel and previous) contains a change (the
>>     line "long ran_x[KK];" is commented, and ran_x is defined elsewhere
>>     in the program). The patch contains the _same_ change; I let you
>>     decide if this is acceptable in regard to Knuth's request.

--- RNG.c.old	Fri Mar  1 08:42:45 2002
+++ RNG.c	Fri Mar  1 14:26:08 2002
@@ -543,6 +543,7 @@
 #define ran_x             dummy
 
 /* ===================  Knuth TAOCP ========================== */
+/* From http://www-cs-faculty.Stanford.EDU/~knuth/programs/rng.c */
 
          /* Please do NOT alter this part of the code */
 
@@ -551,15 +552,17 @@
  *    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 pages 171 and following of Volume 2).              */
+ *    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. */
 
 /*    If you find any bugs, please report them immediately to
  *                 taocp@cs.stanford.edu
  *    (and you will be rewarded if the bug is genuine). Thanks!            */
 
 /************ see the book for explanations and caveats! *******************/
-
-/* the old C calling conventions are used here, for reasons of portability */
+/************ in particular, you need two's complement arithmetic **********/
 
 #define KK 100                     /* the long lag */
 #define LL  37                     /* the short lag */
@@ -568,10 +571,13 @@
 
 /*long ran_x[KK]; */                   /* the generator state */
 
-/* void ran_array(long aa[],int n) */
+#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];
@@ -580,57 +586,57 @@
   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_sentinel=-1;
+long *ran_arr_ptr=&ran_arr_sentinel; /* the next random number, or -1 */
+
+#define ran_arr_next() (*ran_arr_ptr>=0? *ran_arr_ptr++: ran_arr_cycle())
+long ran_arr_cycle()
+{
+  ran_array(ran_arr_buf,QUALITY);
+  ran_arr_buf[100]=-1;
+  ran_arr_ptr=ran_arr_buf+1;
+  return ran_arr_buf[0];
+}
+
 #define TT  70   /* guaranteed separation between streams */
 #define is_odd(x)  ((x)&1)          /* units bit of x */
-#define evenize(x) ((x)&(MM-2))   /* make x even */
 
-/* void ran_start(long seed) */
+#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=evenize(seed+2);
+  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 */
   }
-  for (;j<KK+KK-1;j++) x[j]=0;
   x[1]++;              /* make x[1] (and only x[1]) odd */
-  ss=seed&(MM-1);
-  t=TT-1; while (t) {
-    for (j=KK-1;j>0;j--) x[j+j]=x[j];  /* "square" */
-    for (j=KK+KK-2;j>KK-LL;j-=2) x[KK+KK-1-j]=evenize(x[j]);
-    for (j=KK+KK-2;j>=KK;j--) if(is_odd(x[j])) {
-      x[j-(KK-LL)]=mod_diff(x[j-(KK-LL)],x[j]);
+  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 */
-      if (is_odd(x[KK])) x[LL]=mod_diff(x[LL],x[KK]);
+      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];
-}
-
-/* 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_sentinel=-1;
-long *ran_arr_ptr=&ran_arr_sentinel; /* the next random number, or -1 */
-
-#define ran_arr_next() (*ran_arr_ptr>=0? *ran_arr_ptr++: ran_arr_cycle())
-long ran_arr_cycle()
-{
-  ran_array(ran_arr_buf,QUALITY);
-  ran_arr_buf[100]=-1;
-  ran_arr_ptr=ran_arr_buf+1;
-  return ran_arr_buf[0];
+  for (j=0;j<10;j++) ran_array(x,KK+KK-1); /* warm things up */
+  ran_arr_ptr=&ran_arr_sentinel;
 }
 
 /* ===================== end of Knuth's code ====================== */


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._


[Rd] Weakness in Knuth-TAOCP RNG (fwd) (PR#1336)
schutz@wehi.edu.au