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

@@ -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

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

noao/onedspec/splot/spdeblend.x

@@ -1,3 +1,4 @@
+include	<math.h>
 include	<error.h>
 include	<gset.h>
 
@@ -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

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