GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
mltmod.f
Go to the documentation of this file.
1  INTEGER FUNCTION mltmod(a,s,m)
2 C**********************************************************************
3 C
4 C INTEGER 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 h
25  parameter(h=32768)
26 C ..
27 C .. Scalar Arguments ..
28  INTEGER a,m,s
29 C ..
30 C .. Local Scalars ..
31  INTEGER 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