GNU Octave
4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
Main Page
Namespaces
Classes
Files
File List
File Members
All
Classes
Namespaces
Files
Functions
Variables
Typedefs
Enumerations
Enumerator
Properties
Friends
Macros
Pages
liboctave
cruft
ranlib
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
to
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
Definition:
Quad-opts.cc:233
Generated on Wed May 10 2017 15:42:57 for GNU Octave by
1.8.8