thr3ads.net - R devel - [Rd] Weakness in Knuth-TAOCP RNG (fwd) (PR#1336) [Mar 2002]

If this information is useful, please help other people find it:
Share via:

schutz@wehi.edu.au

2002-Mar-01 06:06 UTC

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

This message is in MIME format.  The first part should be readable text,
  while the remaining parts are likely unreadable without MIME-aware tools.
  Send mail to mime@docserver.cac.washington.edu for more info.

---559023410-1141662977-1014960253=:3119
Content-Type: TEXT/PLAIN; CHARSET=iso-8859-1
Content-Transfer-Encoding: QUOTED-PRINTABLE
Content-ID: <Pine.SOL.4.33.0203011649142.3119@unix09.wehi.edu.au>

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.

Frederic
--
Fr=E9d=E9ric SCH=DCTZ                      email: schutz@wehi.edu.au
Genetics and Bioinformatics Group
The Walter and Eliza Hall Institute of Medical Research
Phone: +61 3 9345 2627               Fax: +61 3 9347 0852
Post Office
Royal Melbourne Hospital, 3050
Victoria, Australia

---559023410-1141662977-1014960253=:3119
Content-Type: TEXT/PLAIN; CHARSET=US-ASCII; NAME="RNG.c.patch"
Content-Transfer-Encoding: BASE64
Content-ID: <Pine.SOL.4.33.0203011624131.3119@unix09.wehi.edu.au>
Content-Description: 
Content-Disposition: ATTACHMENT; FILENAME="RNG.c.patch"

LS0tIFJORy5jLm9sZAlGcmkgTWFyICAxIDA4OjQyOjQ1IDIwMDINCisrKyBS
TkcuYwlGcmkgTWFyICAxIDE0OjI2OjA4IDIwMDINCkBAIC01NDMsNiArNTQz
LDcgQEANCiAjZGVmaW5lIHJhbl94ICAgICAgICAgICAgIGR1bW15DQogDQog
LyogPT09PT09PT09PT09PT09PT09PSAgS251dGggVEFPQ1AgPT09PT09PT09
PT09PT09PT09PT09PT09PT0gKi8NCisvKiBGcm9tIGh0dHA6Ly93d3ctY3Mt
ZmFjdWx0eS5TdGFuZm9yZC5FRFUvfmtudXRoL3Byb2dyYW1zL3JuZy5jICov
DQogDQogICAgICAgICAgLyogUGxlYXNlIGRvIE5PVCBhbHRlciB0aGlzIHBh
cnQgb2YgdGhlIGNvZGUgKi8NCiANCkBAIC01NTEsMTUgKzU1MiwxNyBAQA0K
ICAqICAgIEl0IGlzIGV4cGxhaW5lZCBpbiBTZW1pbnVtZXJpY2FsIEFsZ29y
aXRobXMsIDNyZCBlZGl0aW9uLCBTZWN0aW9uIDMuNg0KICAqICAgIChvciBp
biB0aGUgZXJyYXRhIHRvIHRoZSAybmQgZWRpdGlvbiAtLS0gc2VlDQogICog
ICAgICAgIGh0dHA6Ly93d3ctY3MtZmFjdWx0eS5zdGFuZm9yZC5lZHUvfmtu
dXRoL3Rhb2NwLmh0bWwNCi0gKiAgICBpbiB0aGUgY2hhbmdlcyB0byBwYWdl
cyAxNzEgYW5kIGZvbGxvd2luZyBvZiBWb2x1bWUgMikuICAgICAgICAgICAg
ICAqLw0KKyAqICAgIGluIHRoZSBjaGFuZ2VzIHRvIFZvbHVtZSAyIG9uIHBh
Z2VzIDE3MSBhbmQgZm9sbG93aW5nKS4gICAgICAgICAgICAgICovDQorDQor
LyogICAgTi5CLiBUaGUgTU9ESUZJQ0FUSU9OUyBpbnRyb2R1Y2VkIGluIHRo
ZSA5dGggcHJpbnRpbmcgKDIwMDIpIGFyZQ0KKyAgICAgIGluY2x1ZGVkIGhl
cmU7IHRoZXJlJ3Mgbm8gYmFja3dhcmRzIGNvbXBhdGliaWxpdHkgd2l0aCB0
aGUgb3JpZ2luYWwuICovDQogDQogLyogICAgSWYgeW91IGZpbmQgYW55IGJ1
Z3MsIHBsZWFzZSByZXBvcnQgdGhlbSBpbW1lZGlhdGVseSB0bw0KICAqICAg
ICAgICAgICAgICAgICB0YW9jcEBjcy5zdGFuZm9yZC5lZHUNCiAgKiAgICAo
YW5kIHlvdSB3aWxsIGJlIHJld2FyZGVkIGlmIHRoZSBidWcgaXMgZ2VudWlu
ZSkuIFRoYW5rcyEgICAgICAgICAgICAqLw0KIA0KIC8qKioqKioqKioqKiog
c2VlIHRoZSBib29rIGZvciBleHBsYW5hdGlvbnMgYW5kIGNhdmVhdHMhICoq
KioqKioqKioqKioqKioqKiovDQotDQotLyogdGhlIG9sZCBDIGNhbGxpbmcg
Y29udmVudGlvbnMgYXJlIHVzZWQgaGVyZSwgZm9yIHJlYXNvbnMgb2YgcG9y
dGFiaWxpdHkgKi8NCisvKioqKioqKioqKioqIGluIHBhcnRpY3VsYXIsIHlv
dSBuZWVkIHR3bydzIGNvbXBsZW1lbnQgYXJpdGhtZXRpYyAqKioqKioqKioq
Lw0KIA0KICNkZWZpbmUgS0sgMTAwICAgICAgICAgICAgICAgICAgICAgLyog
dGhlIGxvbmcgbGFnICovDQogI2RlZmluZSBMTCAgMzcgICAgICAgICAgICAg
ICAgICAgICAvKiB0aGUgc2hvcnQgbGFnICovDQpAQCAtNTY4LDEwICs1NzEs
MTMgQEANCiANCiAvKmxvbmcgcmFuX3hbS0tdOyAqLyAgICAgICAgICAgICAg
ICAgICAvKiB0aGUgZ2VuZXJhdG9yIHN0YXRlICovDQogDQotLyogdm9pZCBy
YW5fYXJyYXkobG9uZyBhYVtdLGludCBuKSAqLw0KKyNpZmRlZiBfX1NURENf
Xw0KK3ZvaWQgcmFuX2FycmF5KGxvbmcgYWFbXSxpbnQgbikNCisjZWxzZQ0K
IHZvaWQgcmFuX2FycmF5KGFhLG4pICAgIC8qIHB1dCBuIG5ldyByYW5kb20g
bnVtYmVycyBpbiBhYSAqLw0KICAgbG9uZyAqYWE7ICAgLyogZGVzdGluYXRp
b24gKi8NCiAgIGludCBuOyAgICAgIC8qIGFycmF5IGxlbmd0aCAobXVzdCBi
ZSBhdCBsZWFzdCBLSykgKi8NCisjZW5kaWYNCiB7DQogICByZWdpc3RlciBp
bnQgaSxqOw0KICAgZm9yIChqPTA7ajxLSztqKyspIGFhW2pdPXJhbl94W2pd
Ow0KQEAgLTU4MCw1NyArNTg2LDU3IEBADQogICBmb3IgKDtpPEtLO2krKyxq
KyspIHJhbl94W2ldPW1vZF9kaWZmKGFhW2otS0tdLHJhbl94W2ktTExdKTsN
CiB9DQogDQorLyogdGhlIGZvbGxvd2luZyByb3V0aW5lcyBhcmUgZnJvbSBl
eGVyY2lzZSAzLjYtLTE1ICovDQorLyogYWZ0ZXIgY2FsbGluZyByYW5fc3Rh
cnQsIGdldCBuZXcgcmFuZG9tcyBieSwgZS5nLiwgIng9cmFuX2Fycl9uZXh0
KCkiICovDQorDQorI2RlZmluZSBRVUFMSVRZIDEwMDkgLyogcmVjb21tZW5k
ZWQgcXVhbGl0eSBsZXZlbCBmb3IgaGlnaC1yZXMgdXNlICovDQorbG9uZyBy
YW5fYXJyX2J1ZltRVUFMSVRZXTsNCitsb25nIHJhbl9hcnJfc2VudGluZWw9
LTE7DQorbG9uZyAqcmFuX2Fycl9wdHI9JnJhbl9hcnJfc2VudGluZWw7IC8q
IHRoZSBuZXh0IHJhbmRvbSBudW1iZXIsIG9yIC0xICovDQorDQorI2RlZmlu
ZSByYW5fYXJyX25leHQoKSAoKnJhbl9hcnJfcHRyPj0wPyAqcmFuX2Fycl9w
dHIrKzogcmFuX2Fycl9jeWNsZSgpKQ0KK2xvbmcgcmFuX2Fycl9jeWNsZSgp
DQorew0KKyAgcmFuX2FycmF5KHJhbl9hcnJfYnVmLFFVQUxJVFkpOw0KKyAg
cmFuX2Fycl9idWZbMTAwXT0tMTsNCisgIHJhbl9hcnJfcHRyPXJhbl9hcnJf
YnVmKzE7DQorICByZXR1cm4gcmFuX2Fycl9idWZbMF07DQorfQ0KKw0KICNk
ZWZpbmUgVFQgIDcwICAgLyogZ3VhcmFudGVlZCBzZXBhcmF0aW9uIGJldHdl
ZW4gc3RyZWFtcyAqLw0KICNkZWZpbmUgaXNfb2RkKHgpICAoKHgpJjEpICAg
ICAgICAgIC8qIHVuaXRzIGJpdCBvZiB4ICovDQotI2RlZmluZSBldmVuaXpl
KHgpICgoeCkmKE1NLTIpKSAgIC8qIG1ha2UgeCBldmVuICovDQogDQotLyog
dm9pZCByYW5fc3RhcnQobG9uZyBzZWVkKSAqLw0KKyNpZmRlZiBfX1NURENf
Xw0KK3ZvaWQgcmFuX3N0YXJ0KGxvbmcgc2VlZCkNCisjZWxzZQ0KIHZvaWQg
cmFuX3N0YXJ0KHNlZWQpICAgIC8qIGRvIHRoaXMgYmVmb3JlIHVzaW5nIHJh
bl9hcnJheSAqLw0KICAgbG9uZyBzZWVkOyAgICAgICAgICAgIC8qIHNlbGVj
dG9yIGZvciBkaWZmZXJlbnQgc3RyZWFtcyAqLw0KKyNlbmRpZg0KIHsNCiAg
IHJlZ2lzdGVyIGludCB0LGo7DQogICBsb25nIHhbS0srS0stMV07ICAgICAg
ICAgICAgICAvKiB0aGUgcHJlcGFyYXRpb24gYnVmZmVyICovDQotICByZWdp
c3RlciBsb25nIHNzPWV2ZW5pemUoc2VlZCsyKTsNCisgIHJlZ2lzdGVyIGxv
bmcgc3M9KHNlZWQrMikmKE1NLTIpOw0KICAgZm9yIChqPTA7ajxLSztqKysp
IHsNCiAgICAgeFtqXT1zczsgICAgICAgICAgICAgICAgICAgICAgLyogYm9v
dHN0cmFwIHRoZSBidWZmZXIgKi8NCiAgICAgc3M8PD0xOyBpZiAoc3M+PU1N
KSBzcy09TU0tMjsgLyogY3ljbGljIHNoaWZ0IDI5IGJpdHMgKi8NCiAgIH0N
Ci0gIGZvciAoO2o8S0srS0stMTtqKyspIHhbal09MDsNCiAgIHhbMV0rKzsg
ICAgICAgICAgICAgIC8qIG1ha2UgeFsxXSAoYW5kIG9ubHkgeFsxXSkgb2Rk
ICovDQotICBzcz1zZWVkJihNTS0xKTsNCi0gIHQ9VFQtMTsgd2hpbGUgKHQp
IHsNCi0gICAgZm9yIChqPUtLLTE7aj4wO2otLSkgeFtqK2pdPXhbal07ICAv
KiAic3F1YXJlIiAqLw0KLSAgICBmb3IgKGo9S0srS0stMjtqPktLLUxMO2ot
PTIpIHhbS0srS0stMS1qXT1ldmVuaXplKHhbal0pOw0KLSAgICBmb3IgKGo9
S0srS0stMjtqPj1LSztqLS0pIGlmKGlzX29kZCh4W2pdKSkgew0KLSAgICAg
IHhbai0oS0stTEwpXT1tb2RfZGlmZih4W2otKEtLLUxMKV0seFtqXSk7DQor
ICBmb3IgKHNzPXNlZWQmKE1NLTEpLHQ9VFQtMTsgdDsgKSB7ICAgICAgIA0K
KyAgICBmb3IgKGo9S0stMTtqPjA7ai0tKSB4W2oral09eFtqXSwgeFtqK2ot
MV09MDsgLyogInNxdWFyZSIgKi8NCisgICAgZm9yIChqPUtLK0tLLTI7aj49
S0s7ai0tKQ0KKyAgICAgIHhbai0oS0stTEwpXT1tb2RfZGlmZih4W2otKEtL
LUxMKV0seFtqXSksDQogICAgICAgeFtqLUtLXT1tb2RfZGlmZih4W2otS0td
LHhbal0pOw0KLSAgICB9DQogICAgIGlmIChpc19vZGQoc3MpKSB7ICAgICAg
ICAgICAgICAvKiAibXVsdGlwbHkgYnkgeiIgKi8NCiAgICAgICBmb3IgKGo9
S0s7aj4wO2otLSkgIHhbal09eFtqLTFdOw0KICAgICAgIHhbMF09eFtLS107
ICAgICAgICAgICAgLyogc2hpZnQgdGhlIGJ1ZmZlciBjeWNsaWNhbGx5ICov
DQotICAgICAgaWYgKGlzX29kZCh4W0tLXSkpIHhbTExdPW1vZF9kaWZmKHhb
TExdLHhbS0tdKTsNCisgICAgICB4W0xMXT1tb2RfZGlmZih4W0xMXSx4W0tL
XSk7DQogICAgIH0NCiAgICAgaWYgKHNzKSBzcz4+PTE7IGVsc2UgdC0tOw0K
ICAgfQ0KICAgZm9yIChqPTA7ajxMTDtqKyspIHJhbl94W2orS0stTExdPXhb
al07DQogICBmb3IgKDtqPEtLO2orKykgcmFuX3hbai1MTF09eFtqXTsNCi19
DQotDQotLyogdGhlIGZvbGxvd2luZyByb3V0aW5lcyBhcmUgZnJvbSBleGVy
Y2lzZSAzLjYtLTE1ICovDQotLyogYWZ0ZXIgY2FsbGluZyByYW5fc3RhcnQs
IGdldCBuZXcgcmFuZG9tcyBieSwgZS5nLiwgIng9cmFuX2Fycl9uZXh0KCki
ICovDQotDQotI2RlZmluZSBRVUFMSVRZIDEwMDkgLyogcmVjb21tZW5kZWQg
cXVhbGl0eSBsZXZlbCBmb3IgaGlnaC1yZXMgdXNlICovDQotbG9uZyByYW5f
YXJyX2J1ZltRVUFMSVRZXTsNCi1sb25nIHJhbl9hcnJfc2VudGluZWw9LTE7
DQotbG9uZyAqcmFuX2Fycl9wdHI9JnJhbl9hcnJfc2VudGluZWw7IC8qIHRo
ZSBuZXh0IHJhbmRvbSBudW1iZXIsIG9yIC0xICovDQotDQotI2RlZmluZSBy
YW5fYXJyX25leHQoKSAoKnJhbl9hcnJfcHRyPj0wPyAqcmFuX2Fycl9wdHIr
KzogcmFuX2Fycl9jeWNsZSgpKQ0KLWxvbmcgcmFuX2Fycl9jeWNsZSgpDQot
ew0KLSAgcmFuX2FycmF5KHJhbl9hcnJfYnVmLFFVQUxJVFkpOw0KLSAgcmFu
X2Fycl9idWZbMTAwXT0tMTsNCi0gIHJhbl9hcnJfcHRyPXJhbl9hcnJfYnVm
KzE7DQotICByZXR1cm4gcmFuX2Fycl9idWZbMF07DQorICBmb3IgKGo9MDtq
PDEwO2orKykgcmFuX2FycmF5KHgsS0srS0stMSk7IC8qIHdhcm0gdGhpbmdz
IHVwICovDQorICByYW5fYXJyX3B0cj0mcmFuX2Fycl9zZW50aW5lbDsNCiB9
DQogDQogLyogPT09PT09PT09PT09PT09PT09PT09IGVuZCBvZiBLbnV0aCdz
IGNvZGUgPT09PT09PT09PT09PT09PT09PT09PSAqLw0K
---559023410-1141662977-1014960253=:3119--

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

ripley@stats.ox.ac.uk

2002-Mar-01 06:38 UTC

head link

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

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?

On Fri, 1 Mar 2002 schutz@wehi.edu.au wrote:
>   This message is in MIME format.  The first part should be readable text,
>   while the remaining parts are likely unreadable without MIME-aware tools.
>   Send mail to mime@docserver.cac.washington.edu for more info.
>
> ---559023410-1141662977-1014960253=:3119
> Content-Type: TEXT/PLAIN; CHARSET=iso-8859-1
> Content-Transfer-Encoding: QUOTED-PRINTABLE
> Content-ID: <Pine.SOL.4.33.0203011649142.3119@unix09.wehi.edu.au>
>
> 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.
>
> Frederic
> --
> Fr=E9d=E9ric SCH=DCTZ                      email: schutz@wehi.edu.au
> Genetics and Bioinformatics Group
> The Walter and Eliza Hall Institute of Medical Research
> Phone: +61 3 9345 2627               Fax: +61 3 9347 0852
> Post Office
> Royal Melbourne Hospital, 3050
> Victoria, Australia
>
> ---559023410-1141662977-1014960253=:3119
> Content-Type: TEXT/PLAIN; CHARSET=US-ASCII; NAME="RNG.c.patch"
> Content-Transfer-Encoding: BASE64
> Content-ID: <Pine.SOL.4.33.0203011624131.3119@unix09.wehi.edu.au>
> Content-Description:
> Content-Disposition: ATTACHMENT; FILENAME="RNG.c.patch"
>
> LS0tIFJORy5jLm9sZAlGcmkgTWFyICAxIDA4OjQyOjQ1IDIwMDINCisrKyBS
> TkcuYwlGcmkgTWFyICAxIDE0OjI2OjA4IDIwMDINCkBAIC01NDMsNiArNTQz
> LDcgQEANCiAjZGVmaW5lIHJhbl94ICAgICAgICAgICAgIGR1bW15DQogDQog
> LyogPT09PT09PT09PT09PT09PT09PSAgS251dGggVEFPQ1AgPT09PT09PT09
> PT09PT09PT09PT09PT09PT0gKi8NCisvKiBGcm9tIGh0dHA6Ly93d3ctY3Mt
> ZmFjdWx0eS5TdGFuZm9yZC5FRFUvfmtudXRoL3Byb2dyYW1zL3JuZy5jICov
> DQogDQogICAgICAgICAgLyogUGxlYXNlIGRvIE5PVCBhbHRlciB0aGlzIHBh
> cnQgb2YgdGhlIGNvZGUgKi8NCiANCkBAIC01NTEsMTUgKzU1MiwxNyBAQA0K
> ICAqICAgIEl0IGlzIGV4cGxhaW5lZCBpbiBTZW1pbnVtZXJpY2FsIEFsZ29y
> aXRobXMsIDNyZCBlZGl0aW9uLCBTZWN0aW9uIDMuNg0KICAqICAgIChvciBp
> biB0aGUgZXJyYXRhIHRvIHRoZSAybmQgZWRpdGlvbiAtLS0gc2VlDQogICog
> ICAgICAgIGh0dHA6Ly93d3ctY3MtZmFjdWx0eS5zdGFuZm9yZC5lZHUvfmtu
> dXRoL3Rhb2NwLmh0bWwNCi0gKiAgICBpbiB0aGUgY2hhbmdlcyB0byBwYWdl
> cyAxNzEgYW5kIGZvbGxvd2luZyBvZiBWb2x1bWUgMikuICAgICAgICAgICAg
> ICAqLw0KKyAqICAgIGluIHRoZSBjaGFuZ2VzIHRvIFZvbHVtZSAyIG9uIHBh
> Z2VzIDE3MSBhbmQgZm9sbG93aW5nKS4gICAgICAgICAgICAgICovDQorDQor
> LyogICAgTi5CLiBUaGUgTU9ESUZJQ0FUSU9OUyBpbnRyb2R1Y2VkIGluIHRo
> ZSA5dGggcHJpbnRpbmcgKDIwMDIpIGFyZQ0KKyAgICAgIGluY2x1ZGVkIGhl
> cmU7IHRoZXJlJ3Mgbm8gYmFja3dhcmRzIGNvbXBhdGliaWxpdHkgd2l0aCB0
> aGUgb3JpZ2luYWwuICovDQogDQogLyogICAgSWYgeW91IGZpbmQgYW55IGJ1
> Z3MsIHBsZWFzZSByZXBvcnQgdGhlbSBpbW1lZGlhdGVseSB0bw0KICAqICAg
> ICAgICAgICAgICAgICB0YW9jcEBjcy5zdGFuZm9yZC5lZHUNCiAgKiAgICAo
> YW5kIHlvdSB3aWxsIGJlIHJld2FyZGVkIGlmIHRoZSBidWcgaXMgZ2VudWlu
> ZSkuIFRoYW5rcyEgICAgICAgICAgICAqLw0KIA0KIC8qKioqKioqKioqKiog
> c2VlIHRoZSBib29rIGZvciBleHBsYW5hdGlvbnMgYW5kIGNhdmVhdHMhICoq
> KioqKioqKioqKioqKioqKiovDQotDQotLyogdGhlIG9sZCBDIGNhbGxpbmcg
> Y29udmVudGlvbnMgYXJlIHVzZWQgaGVyZSwgZm9yIHJlYXNvbnMgb2YgcG9y
> dGFiaWxpdHkgKi8NCisvKioqKioqKioqKioqIGluIHBhcnRpY3VsYXIsIHlv
> dSBuZWVkIHR3bydzIGNvbXBsZW1lbnQgYXJpdGhtZXRpYyAqKioqKioqKioq
> Lw0KIA0KICNkZWZpbmUgS0sgMTAwICAgICAgICAgICAgICAgICAgICAgLyog
> dGhlIGxvbmcgbGFnICovDQogI2RlZmluZSBMTCAgMzcgICAgICAgICAgICAg
> ICAgICAgICAvKiB0aGUgc2hvcnQgbGFnICovDQpAQCAtNTY4LDEwICs1NzEs
> MTMgQEANCiANCiAvKmxvbmcgcmFuX3hbS0tdOyAqLyAgICAgICAgICAgICAg
> ICAgICAvKiB0aGUgZ2VuZXJhdG9yIHN0YXRlICovDQogDQotLyogdm9pZCBy
> YW5fYXJyYXkobG9uZyBhYVtdLGludCBuKSAqLw0KKyNpZmRlZiBfX1NURENf
> Xw0KK3ZvaWQgcmFuX2FycmF5KGxvbmcgYWFbXSxpbnQgbikNCisjZWxzZQ0K
> IHZvaWQgcmFuX2FycmF5KGFhLG4pICAgIC8qIHB1dCBuIG5ldyByYW5kb20g
> bnVtYmVycyBpbiBhYSAqLw0KICAgbG9uZyAqYWE7ICAgLyogZGVzdGluYXRp
> b24gKi8NCiAgIGludCBuOyAgICAgIC8qIGFycmF5IGxlbmd0aCAobXVzdCBi
> ZSBhdCBsZWFzdCBLSykgKi8NCisjZW5kaWYNCiB7DQogICByZWdpc3RlciBp
> bnQgaSxqOw0KICAgZm9yIChqPTA7ajxLSztqKyspIGFhW2pdPXJhbl94W2pd
> Ow0KQEAgLTU4MCw1NyArNTg2LDU3IEBADQogICBmb3IgKDtpPEtLO2krKyxq
> KyspIHJhbl94W2ldPW1vZF9kaWZmKGFhW2otS0tdLHJhbl94W2ktTExdKTsN
> CiB9DQogDQorLyogdGhlIGZvbGxvd2luZyByb3V0aW5lcyBhcmUgZnJvbSBl
> eGVyY2lzZSAzLjYtLTE1ICovDQorLyogYWZ0ZXIgY2FsbGluZyByYW5fc3Rh
> cnQsIGdldCBuZXcgcmFuZG9tcyBieSwgZS5nLiwgIng9cmFuX2Fycl9uZXh0
> KCkiICovDQorDQorI2RlZmluZSBRVUFMSVRZIDEwMDkgLyogcmVjb21tZW5k
> ZWQgcXVhbGl0eSBsZXZlbCBmb3IgaGlnaC1yZXMgdXNlICovDQorbG9uZyBy
> YW5fYXJyX2J1ZltRVUFMSVRZXTsNCitsb25nIHJhbl9hcnJfc2VudGluZWw9
> LTE7DQorbG9uZyAqcmFuX2Fycl9wdHI9JnJhbl9hcnJfc2VudGluZWw7IC8q
> IHRoZSBuZXh0IHJhbmRvbSBudW1iZXIsIG9yIC0xICovDQorDQorI2RlZmlu
> ZSByYW5fYXJyX25leHQoKSAoKnJhbl9hcnJfcHRyPj0wPyAqcmFuX2Fycl9w
> dHIrKzogcmFuX2Fycl9jeWNsZSgpKQ0KK2xvbmcgcmFuX2Fycl9jeWNsZSgp
> DQorew0KKyAgcmFuX2FycmF5KHJhbl9hcnJfYnVmLFFVQUxJVFkpOw0KKyAg
> cmFuX2Fycl9idWZbMTAwXT0tMTsNCisgIHJhbl9hcnJfcHRyPXJhbl9hcnJf
> YnVmKzE7DQorICByZXR1cm4gcmFuX2Fycl9idWZbMF07DQorfQ0KKw0KICNk
> ZWZpbmUgVFQgIDcwICAgLyogZ3VhcmFudGVlZCBzZXBhcmF0aW9uIGJldHdl
> ZW4gc3RyZWFtcyAqLw0KICNkZWZpbmUgaXNfb2RkKHgpICAoKHgpJjEpICAg
> ICAgICAgIC8qIHVuaXRzIGJpdCBvZiB4ICovDQotI2RlZmluZSBldmVuaXpl
> KHgpICgoeCkmKE1NLTIpKSAgIC8qIG1ha2UgeCBldmVuICovDQogDQotLyog
> dm9pZCByYW5fc3RhcnQobG9uZyBzZWVkKSAqLw0KKyNpZmRlZiBfX1NURENf
> Xw0KK3ZvaWQgcmFuX3N0YXJ0KGxvbmcgc2VlZCkNCisjZWxzZQ0KIHZvaWQg
> cmFuX3N0YXJ0KHNlZWQpICAgIC8qIGRvIHRoaXMgYmVmb3JlIHVzaW5nIHJh
> bl9hcnJheSAqLw0KICAgbG9uZyBzZWVkOyAgICAgICAgICAgIC8qIHNlbGVj
> dG9yIGZvciBkaWZmZXJlbnQgc3RyZWFtcyAqLw0KKyNlbmRpZg0KIHsNCiAg
> IHJlZ2lzdGVyIGludCB0LGo7DQogICBsb25nIHhbS0srS0stMV07ICAgICAg
> ICAgICAgICAvKiB0aGUgcHJlcGFyYXRpb24gYnVmZmVyICovDQotICByZWdp
> c3RlciBsb25nIHNzPWV2ZW5pemUoc2VlZCsyKTsNCisgIHJlZ2lzdGVyIGxv
> bmcgc3M9KHNlZWQrMikmKE1NLTIpOw0KICAgZm9yIChqPTA7ajxLSztqKysp
> IHsNCiAgICAgeFtqXT1zczsgICAgICAgICAgICAgICAgICAgICAgLyogYm9v
> dHN0cmFwIHRoZSBidWZmZXIgKi8NCiAgICAgc3M8PD0xOyBpZiAoc3M+PU1N
> KSBzcy09TU0tMjsgLyogY3ljbGljIHNoaWZ0IDI5IGJpdHMgKi8NCiAgIH0N
> Ci0gIGZvciAoO2o8S0srS0stMTtqKyspIHhbal09MDsNCiAgIHhbMV0rKzsg
> ICAgICAgICAgICAgIC8qIG1ha2UgeFsxXSAoYW5kIG9ubHkgeFsxXSkgb2Rk
> ICovDQotICBzcz1zZWVkJihNTS0xKTsNCi0gIHQ9VFQtMTsgd2hpbGUgKHQp
> IHsNCi0gICAgZm9yIChqPUtLLTE7aj4wO2otLSkgeFtqK2pdPXhbal07ICAv
> KiAic3F1YXJlIiAqLw0KLSAgICBmb3IgKGo9S0srS0stMjtqPktLLUxMO2ot
> PTIpIHhbS0srS0stMS1qXT1ldmVuaXplKHhbal0pOw0KLSAgICBmb3IgKGo9
> S0srS0stMjtqPj1LSztqLS0pIGlmKGlzX29kZCh4W2pdKSkgew0KLSAgICAg
> IHhbai0oS0stTEwpXT1tb2RfZGlmZih4W2otKEtLLUxMKV0seFtqXSk7DQor
> ICBmb3IgKHNzPXNlZWQmKE1NLTEpLHQ9VFQtMTsgdDsgKSB7ICAgICAgIA0K
> KyAgICBmb3IgKGo9S0stMTtqPjA7ai0tKSB4W2oral09eFtqXSwgeFtqK2ot
> MV09MDsgLyogInNxdWFyZSIgKi8NCisgICAgZm9yIChqPUtLK0tLLTI7aj49
> S0s7ai0tKQ0KKyAgICAgIHhbai0oS0stTEwpXT1tb2RfZGlmZih4W2otKEtL
> LUxMKV0seFtqXSksDQogICAgICAgeFtqLUtLXT1tb2RfZGlmZih4W2otS0td
> LHhbal0pOw0KLSAgICB9DQogICAgIGlmIChpc19vZGQoc3MpKSB7ICAgICAg
> ICAgICAgICAvKiAibXVsdGlwbHkgYnkgeiIgKi8NCiAgICAgICBmb3IgKGo9
> S0s7aj4wO2otLSkgIHhbal09eFtqLTFdOw0KICAgICAgIHhbMF09eFtLS107
> ICAgICAgICAgICAgLyogc2hpZnQgdGhlIGJ1ZmZlciBjeWNsaWNhbGx5ICov
> DQotICAgICAgaWYgKGlzX29kZCh4W0tLXSkpIHhbTExdPW1vZF9kaWZmKHhb
> TExdLHhbS0tdKTsNCisgICAgICB4W0xMXT1tb2RfZGlmZih4W0xMXSx4W0tL
> XSk7DQogICAgIH0NCiAgICAgaWYgKHNzKSBzcz4+PTE7IGVsc2UgdC0tOw0K
> ICAgfQ0KICAgZm9yIChqPTA7ajxMTDtqKyspIHJhbl94W2orS0stTExdPXhb
> al07DQogICBmb3IgKDtqPEtLO2orKykgcmFuX3hbai1MTF09eFtqXTsNCi19
> DQotDQotLyogdGhlIGZvbGxvd2luZyByb3V0aW5lcyBhcmUgZnJvbSBleGVy
> Y2lzZSAzLjYtLTE1ICovDQotLyogYWZ0ZXIgY2FsbGluZyByYW5fc3RhcnQs
> IGdldCBuZXcgcmFuZG9tcyBieSwgZS5nLiwgIng9cmFuX2Fycl9uZXh0KCki
> ICovDQotDQotI2RlZmluZSBRVUFMSVRZIDEwMDkgLyogcmVjb21tZW5kZWQg
> cXVhbGl0eSBsZXZlbCBmb3IgaGlnaC1yZXMgdXNlICovDQotbG9uZyByYW5f
> YXJyX2J1ZltRVUFMSVRZXTsNCi1sb25nIHJhbl9hcnJfc2VudGluZWw9LTE7
> DQotbG9uZyAqcmFuX2Fycl9wdHI9JnJhbl9hcnJfc2VudGluZWw7IC8qIHRo
> ZSBuZXh0IHJhbmRvbSBudW1iZXIsIG9yIC0xICovDQotDQotI2RlZmluZSBy
> YW5fYXJyX25leHQoKSAoKnJhbl9hcnJfcHRyPj0wPyAqcmFuX2Fycl9wdHIr
> KzogcmFuX2Fycl9jeWNsZSgpKQ0KLWxvbmcgcmFuX2Fycl9jeWNsZSgpDQot
> ew0KLSAgcmFuX2FycmF5KHJhbl9hcnJfYnVmLFFVQUxJVFkpOw0KLSAgcmFu
> X2Fycl9idWZbMTAwXT0tMTsNCi0gIHJhbl9hcnJfcHRyPXJhbl9hcnJfYnVm
> KzE7DQotICByZXR1cm4gcmFuX2Fycl9idWZbMF07DQorICBmb3IgKGo9MDtq
> PDEwO2orKykgcmFuX2FycmF5KHgsS0srS0stMSk7IC8qIHdhcm0gdGhpbmdz
> IHVwICovDQorICByYW5fYXJyX3B0cj0mcmFuX2Fycl9zZW50aW5lbDsNCiB9
> DQogDQogLyogPT09PT09PT09PT09PT09PT09PT09IGVuZCBvZiBLbnV0aCdz
> IGNvZGUgPT09PT09PT09PT09PT09PT09PT09PSAqLw0K
> ---559023410-1141662977-1014960253=:3119--
>
>
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> 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
>
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>
-- 
Brian D. Ripley,                  ripley@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595


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

schutz@wehi.edu.au

2002-Mar-03 03:17 UTC

head link

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

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

Maybe Matching Threads

algorithm reference for sample() - Knuth

R devel - Mar 2002 - Weakness in Knuth-TAOCP RNG (fwd) (PR#1336)

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

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

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

Maybe Matching Threads