diff options
Diffstat (limited to 'sys/osb/urand.x')
-rw-r--r-- | sys/osb/urand.x | 55 |
1 files changed, 55 insertions, 0 deletions
diff --git a/sys/osb/urand.x b/sys/osb/urand.x new file mode 100644 index 00000000..84e1fc67 --- /dev/null +++ b/sys/osb/urand.x @@ -0,0 +1,55 @@ +# 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 |