Re: [Rd] Weakness in Knuth-TAOCP RNG (fwd) (PR#1336)

About this list Date view Thread view Subject view Author view Attachment view

From: schutz@wehi.edu.au
Date: Sun 03 Mar 2002 - 14:17:10 EST


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.

Frédéric

---
>> 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._


About this list Date view Thread view Subject view Author view Attachment view

This archive was generated by hypermail 2.1.3 : Wed 16 Oct 2002 - 10:20:50 EST