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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
|
INTEGER FUNCTION mltmod(a,s,m,ierr)
C**********************************************************************
C
C INTEGER FUNCTION MLTMOD(A,S,M)
C
C Returns (A*S) MOD M
C
C This is a transcription from Pascal to Fortran of routine
C MULtMod_Decompos from the paper
C
C L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
C with Splitting Facilities." ACM Transactions on Mathematical
C Software, 17:98-111 (1991)
C
C
C Arguments
C
C
C A, S, M -->
C INTEGER A,S,M
C
C**********************************************************************
include '../stack.h'
C .. Parameters ..
INTEGER h
PARAMETER (h=32768)
C ..
C .. Scalar Arguments ..
INTEGER a,m,s
C ..
C .. Local Scalars ..
INTEGER a0,a1,k,p,q,qh,rh
C ..
C .. Executable Statements ..
C
C H = 2**((b-2)/2) where b = 32 because we are using a 32 bit
C machine. On a different machine recompute H
C
ierr=0
IF (.NOT. (a.LE.0.OR.a.GE.m.OR.s.LE.0.OR.s.GE.m)) GO TO 10
call basout(io,wte,' A, M, S out of order')
call basout(io,wte,'MLTMOD requires: 0 < A < M; 0 < S < M')
mltmod=0
ierr=1
return
10 IF (.NOT. (a.LT.h)) GO TO 20
a0 = a
p = 0
GO TO 120
20 a1 = a/h
a0 = a - h*a1
qh = m/h
rh = m - h*qh
IF (.NOT. (a1.GE.h)) GO TO 50
a1 = a1 - h
k = s/qh
p = h* (s-k*qh) - k*rh
30 IF (.NOT. (p.LT.0)) GO TO 40
p = p + m
GO TO 30
40 GO TO 60
50 p = 0
C
C P = (A2*S*H)MOD M
C
60 IF (.NOT. (a1.NE.0)) GO TO 90
q = m/a1
k = s/q
p = p - k* (m-a1*q)
IF (p.GT.0) p = p - m
p = p + a1* (s-k*q)
70 IF (.NOT. (p.LT.0)) GO TO 80
p = p + m
GO TO 70
80 CONTINUE
90 k = p/qh
C
C P = ((A2*H + A1)*S)MOD M
C
p = h* (p-k*qh) - k*rh
100 IF (.NOT. (p.LT.0)) GO TO 110
p = p + m
GO TO 100
110 CONTINUE
120 IF (.NOT. (a0.NE.0)) GO TO 150
C
C P = ((A2*H + A1)*H*S)MOD M
C
q = m/a0
k = s/q
p = p - k* (m-a0*q)
IF (p.GT.0) p = p - m
p = p + a0* (s-k*q)
130 IF (.NOT. (p.LT.0)) GO TO 140
p = p + m
GO TO 130
140 CONTINUE
150 mltmod = p
C
RETURN
END
|
