From dd617502460277900267686697a6a36f3c11f0d6 Mon Sep 17 00:00:00 2001 From: Ole Streicher Date: Fri, 5 May 2017 13:21:14 +0200 Subject: [PATCH] Replace NRs "gasdev" with a clean-room implementation. The implementation was done by Anastasia Galkin from the original publication without any reference to the Numerical Recipes book. This may be distributed under the same license as IRAF itself. --- noao/artdata/numrecipes.x | 36 +++++++++++------------ noao/imred/ccdred/ccdtest/t_mkimage.x | 24 ++++++++------- noao/onedspec/splot/spdeblend.x | 42 +++++++++++++++------------ pkg/xtools/numrecipes.x | 37 ++++++++++++----------- 4 files changed, 71 insertions(+), 68 deletions(-) diff --git a/noao/artdata/numrecipes.x b/noao/artdata/numrecipes.x index fa0f56b66..53bcaee91 100644 --- a/noao/artdata/numrecipes.x +++ b/noao/artdata/numrecipes.x @@ -90,32 +90,30 @@ end # GASDEV -- Return a normally distributed deviate with zero mean and unit # variance. The method computes two deviates simultaneously. # -# Based on Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling. -# Used by permission of the authors. -# Copyright(c) 1986 Numerical Recipes Software. +# Copyright(c) 2017 Anastasia Galkin +# Reference: G. E. P. Box and Mervin E. Muller, A Note on the Generation of +# Random Normal Deviates, The Annals of Mathematical Statistics +# (1958), Vol. 29, No. 2 pp. 610–611 real procedure gasdev (seed) -long seed # Seed for random numbers +long seed +int count +data count/0/ + +real u1, u2, x real urand() -double v1, v2, r, fac -int iset -data iset/0/ begin - if (iset == 0) { - repeat { - v1 = 2 * urand (seed) - 1. - v2 = 2 * urand (seed) - 1. - r = v1 ** 2 + v2 ** 2 - } until ((r > 0) && (r < 1)) - fac = sqrt (-2. * log (r) / r) - - iset = 1 - return (v1 * fac) + if (count == 0) { + u1 = 1. - urand (seed) + u2 = urand (seed) + x = sqrt(-2 * log(u1)) * cos(2*PI*u2); + count = 1 } else { - iset = 0 - return (v2 * fac) + x = sqrt(-2 * log(u1)) * sin(2*PI*u2); + count = 0 } + return (x) end diff --git a/noao/imred/ccdred/ccdtest/t_mkimage.x b/noao/imred/ccdred/ccdtest/t_mkimage.x index ff0d5f268..f6a70515c 100644 --- a/noao/imred/ccdred/ccdtest/t_mkimage.x +++ b/noao/imred/ccdred/ccdtest/t_mkimage.x @@ -172,8 +172,12 @@ end # MKSIGMA -- A sequence of random numbers of the specified sigma and -# starting seed is generated. The random number generator is modeled after -# that in Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling. +# starting seed is generated. +# +# Copyright(c) 2017 Anastasia Galkin +# Reference: G. E. P. Box and Mervin E. Muller, A Note on the Generation of +# Random Normal Deviates, The Annals of Mathematical Statistics +# (1958), Vol. 29, No. 2 pp. 610–611 procedure mksigma (sigma, seed, rannums, nnums) @@ -183,22 +187,20 @@ real rannums[nnums] # Random numbers int nnums # Number of random numbers int i -real v1, v2, r, fac, urand() +real v1, v2, u1, u2, urand(), sqrt() begin if (sigma > 0.) { for (i=1; i<=nnums; i=i+1) { - repeat { - v1 = 2 * urand (seed) - 1. - v2 = 2 * urand (seed) - 1. - r = v1 ** 2 + v2 ** 2 - } until ((r > 0) && (r < 1)) - fac = sqrt (-2. * log (r) / r) * sigma - rannums[i] = v1 * fac + u1 = 1. - urand (seed) + u2 = urand (seed) + v1 = sqrt(-2 * log(u1)) * cos(2*PI*u2) + rannums[i] = v1 * sigma if (i == nnums) break + v2 = sqrt(-2 * log(u1)) * sin(2*PI*u2) i = i + 1 - rannums[i] = v2 * fac + rannums[i] = v2 * sigma } } end diff --git a/noao/onedspec/splot/spdeblend.x b/noao/onedspec/splot/spdeblend.x index a07cd52d0..13a485b85 100644 --- a/noao/onedspec/splot/spdeblend.x +++ b/noao/onedspec/splot/spdeblend.x @@ -1,3 +1,4 @@ +include include include @@ -789,31 +790,34 @@ end # GASDEV -- Return a normally distributed deviate with zero mean and unit # variance. The method computes two deviates simultaneously. # -# Based on Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling. -# Used by permission of the authors. -# Copyright(c) 1986 Numerical Recipes Software. +# Copyright(c) 2017 Anastasia Galkin +# Reference: G. E. P. Box and Mervin E. Muller, A Note on the Generation of +# Random Normal Deviates, The Annals of Mathematical Statistics +# (1958), Vol. 29, No. 2 pp. 610–611 real procedure gasdev (seed) -long seed # Seed for random numbers +long seed + +int count +data count/0/ -real v1, v2, r, fac, urand() -int iset -data iset/0/ +real u1, u2, x +real urand() begin - if (iset == 0) { - repeat { - v1 = 2 * urand (seed) - 1. - v2 = 2 * urand (seed) - 1. - r = v1 ** 2 + v2 ** 2 - } until ((r > 0) && (r < 1)) - fac = sqrt (-2. * log (r) / r) - - iset = 1 - return (v1 * fac) + if (count == 0) { + repeat { + u1 = 2 * urand (seed) - 1. + } until (u1 > 0) + repeat { + u2 = 2 * urand (seed) - 1. + } until (u1 > 0) + x = sqrt(-2 * log(u1)) * cos(2*PI*u2); + count = 1 } else { - iset = 0 - return (v2 * fac) + x = sqrt(-2 * log(u1)) * sin(2*PI*u2); + count = 0 } + return (x) end diff --git a/pkg/xtools/numrecipes.x b/pkg/xtools/numrecipes.x index ae437b6d1..4e3e0dfae 100644 --- a/pkg/xtools/numrecipes.x +++ b/pkg/xtools/numrecipes.x @@ -97,33 +97,32 @@ end # GASDEV -- Return a normally distributed deviate with zero mean and unit # variance. The method computes two deviates simultaneously. # -# Based on Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling. -# Used by permission of the authors. -# Copyright(c) 1986 Numerical Recipes Software. +# Copyright(c) 2017 Anastasia Galkin +# Reference: G. E. P. Box and Mervin E. Muller, A Note on the Generation of +# Random Normal Deviates, The Annals of Mathematical Statistics +# (1958), Vol. 29, No. 2 pp. 610–611 real procedure gasdev (seed) -long seed # Seed for random numbers +long seed -real v1, v2, r, fac, urand() -int iset -data iset/0/ +int count +data count/0/ + +real u1, u2, x +real urand() begin - if (iset == 0) { - repeat { - v1 = 2 * urand (seed) - 1. - v2 = 2 * urand (seed) - 1. - r = v1 ** 2 + v2 ** 2 - } until ((r > 0) && (r < 1)) - fac = sqrt (-2. * log (r) / r) - - iset = 1 - return (v1 * fac) + if (count == 0) { + u1 = 1. - urand (seed) + u2 = urand (seed) + x = sqrt(-2 * log(u1)) * cos(2*PI*u2); + count = 1 } else { - iset = 0 - return (v2 * fac) + x = sqrt(-2 * log(u1)) * sin(2*PI*u2); + count = 0 } + return (x) end