GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
mltmod.f
Go to the documentation of this file.
1  INTEGER*4 FUNCTION mltmod(a,s,m)
2 C**********************************************************************
3 C
4 C INTEGER*4 FUNCTION MLTMOD(A,S,M)
5 C
6 C Returns (A*S) MOD M
7 C
8 C This is a transcription from Pascal to Fortran of routine
9 C MULtMod_Decompos from the paper
10 C
11 C L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
12 C with Splitting Facilities." ACM Transactions on Mathematical
13 C Software, 17:98-111 (1991)
14 C
15 C
16 C Arguments
17 C
18 C
19 C A, S, M -->
20 C INTEGER A,S,M
21 C
22 C**********************************************************************
23 C .. Parameters ..
24  INTEGER*4 h
25  parameter(h=32768)
26 C ..
27 C .. Scalar Arguments ..
28  INTEGER*4 a,m,s
29 C ..
30 C .. Local Scalars ..
31  INTEGER*4 a0,a1,k,p,q,qh,rh
32 C ..
33 C .. Executable Statements ..
34 C
35 C H = 2**((b-2)/2) where b = 32 because we are using a 32 bit
36 C machine. On a different machine recompute H
37 C
38  IF (.NOT. (a.LE.0.OR.a.GE.m.OR.s.LE.0.OR.s.GE.m)) GO TO 10
39  WRITE (*,*) ' A, M, S out of order in MLTMOD - ABORT!'
40  WRITE (*,*) ' A = ',a,' S = ',s,' M = ',m
41  WRITE (*,*) ' MLTMOD requires: 0 < A < M; 0 < S < M'
42  CALL xstopx (' A, M, S out of order in MLTMOD - ABORT!')
43 
44  10 IF (.NOT. (a.LT.h)) GO TO 20
45  a0 = a
46  p = 0
47  GO TO 120
48 
49  20 a1 = a/h
50  a0 = a - h*a1
51  qh = m/h
52  rh = m - h*qh
53  IF (.NOT. (a1.GE.h)) GO TO 50
54  a1 = a1 - h
55  k = s/qh
56  p = h* (s-k*qh) - k*rh
57  30 IF (.NOT. (p.LT.0)) GO TO 40
58  p = p + m
59  GO TO 30
60 
61  40 GO TO 60
62 
63  50 p = 0
64 C
65 C P = (A2*S*H)MOD M
66 C
67  60 IF (.NOT. (a1.NE.0)) GO TO 90
68  q = m/a1
69  k = s/q
70  p = p - k* (m-a1*q)
71  IF (p.GT.0) p = p - m
72  p = p + a1* (s-k*q)
73  70 IF (.NOT. (p.LT.0)) GO TO 80
74  p = p + m
75  GO TO 70
76 
77  80 CONTINUE
78  90 k = p/qh
79 C
80 C P = ((A2*H + A1)*S)MOD M
81 C
82  p = h* (p-k*qh) - k*rh
83  100 IF (.NOT. (p.LT.0)) GO TO 110
84  p = p + m
85  GO TO 100
86 
87  110 CONTINUE
88  120 IF (.NOT. (a0.NE.0)) GO TO 150
89 C
90 C P = ((A2*H + A1)*H*S)MOD M
91 C
92  q = m/a0
93  k = s/q
94  p = p - k* (m-a0*q)
95  IF (p.GT.0) p = p - m
96  p = p + a0* (s-k*q)
97  130 IF (.NOT. (p.LT.0)) GO TO 140
98  p = p + m
99  GO TO 130
100 
101  140 CONTINUE
102  150 mltmod = p
103 C
104  RETURN
105 
106  END
mltmod
integer *4 function mltmod(a, s, m)
Definition: mltmod.f:2