1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
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
|
