# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.

include	<mach.h>

# URAND -- Universal Random Number Generator.  From "Computer Methods for
# Mathematical Computations", by Forsythe, Malcolm, and Moler, 1977.
# Urand is a uniform random number generator based on theory and suggestions
# given in D.E. Knuth (1969), Vol 2.  Values of URAND will be returned in the
# interval (0,1).  Random numbers are generated by the recursion relation
# (r' = r * a + c) where the art lies in choosing the values for A and C.
#
# [MACHDEP] - NOTE - This routine will not work on machines that do not permit
# integer overflow during multiplication.  In such a case a machine dependent
# routine should be provided in host$as.

real procedure urand (lseed)

long	lseed				# seed value on first call
long	n, a, c, m, mic

real	scale
data	m /0/

int	imul32()

begin
	# When first called, compute multiplier, increment, and miscellaneous
	# constants.

	if (m == 0) {
	    m = MAX_LONG / 2 + 1
	    a = 8 * int (m * atan (1.d0 / 8.d0)) + 5
	    c = 2 * int (m * (0.5d0 - sqrt (3.d0) / 6.d0)) + 1
	    mic = (m - c) + m
	    scale = 0.5 / m
	    lseed = max (1, lseed)
	}

	# Compute next random number, taking care not to cause an arithmetic
	# exception.

	n = imul32 (lseed, a)			# [MACHDEP] - integer overflow
	if (n > mic)
	    n = (n - m) - m
	n = n + c

	if (n / 2 > m)
	    n = (n - m) - m

	if (n < 0)
	    n = (n + m) + m

	lseed = n
	return (n * scale)
end