# SM implementation of algorithm used by Press et al. for deriving
# random numbers with gaussian distribution from random numbers with
# uniform distribution, here the SM built-in function. - bjw, December 2004
# returns a vector of N randoms with rms=1.
# gasdev
gasdev 2 #
set dimen($2) = $1
set _v1 = 2*random($1)-1
set _v2 = 2*random($1)-1
set _rsq = _v1**2 + _v2**2
set _ifgood = (_rsq<1 && _rsq>0) ? 1 : 0
set _fac = (_ifgood==1) ? (-2*ln(_rsq)/_rsq) : 0
set _fac = sqrt(_fac)
set _g1 = _v1 * _fac
set _g2 = _v2 * _fac
# try to fill up with good numbers from the two random vectors
# if the entry wasn't good, we get a number from the second vector.
# because more than half of the numbers fall in the unit circle,
# this should work for any reasonably large N.
define _j 0
do _i = 0, $1-1 {
if ( _ifgood[$_i]==1 ) {
set $2[$_i] = _g1[$_i] } \
else {
while { _ifgood[$_j]==0 } {define _j ($_j+1) }
set $2[$_i] = _g2[$_j]
define _j ($_j+1) }
}
# Copyright 2004, Benjamin Weiner, licensed under GPL.
# You may use or redistribute this code freely, please retain this
# notice. If you incorporate it into a package, you can use the
# package freely for your own research. (You can't sell it without
# distributing the source freely.)