GNU Octave: libcruft/ranlib/mltmod.f Source File

Go to the documentation of this file.
00001       INTEGER FUNCTION mltmod(a,s,m)
00002 C**********************************************************************
00003 C
00004 C     INTEGER FUNCTION MLTMOD(A,S,M)
00005 C
00006 C                    Returns (A*S) MOD M
00007 C
00008 C     This is a transcription from Pascal to Fortran of routine
00009 C     MULtMod_Decompos from the paper
00010 C
00011 C     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
00012 C     with Splitting Facilities." ACM Transactions on Mathematical
00013 C     Software, 17:98-111 (1991)
00014 C
00015 C
00016 C                              Arguments
00017 C
00018 C
00019 C     A, S, M  -->
00020 C                         INTEGER A,S,M
00021 C
00022 C**********************************************************************
00023 C     .. Parameters ..
00024       INTEGER h
00025       PARAMETER (h=32768)
00026 C     ..
00027 C     .. Scalar Arguments ..
00028       INTEGER a,m,s
00029 C     ..
00030 C     .. Local Scalars ..
00031       INTEGER a0,a1,k,p,q,qh,rh
00032 C     ..
00033 C     .. Executable Statements ..
00034 C
00035 C     H = 2**((b-2)/2) where b = 32 because we are using a 32 bit
00036 C      machine. On a different machine recompute H
00037 C
00038       IF (.NOT. (a.LE.0.OR.a.GE.m.OR.s.LE.0.OR.s.GE.m)) GO TO 10
00039       WRITE (*,*) ' A, M, S out of order in MLTMOD - ABORT!'
00040       WRITE (*,*) ' A = ',a,' S = ',s,' M = ',m
00041       WRITE (*,*) ' MLTMOD requires: 0 < A < M; 0 < S < M'
00042       CALL XSTOPX (' A, M, S out of order in MLTMOD - ABORT!')
00043 
00044    10 IF (.NOT. (a.LT.h)) GO TO 20
00045       a0 = a
00046       p = 0
00047       GO TO 120
00048 
00049    20 a1 = a/h
00050       a0 = a - h*a1
00051       qh = m/h
00052       rh = m - h*qh
00053       IF (.NOT. (a1.GE.h)) GO TO 50
00054       a1 = a1 - h
00055       k = s/qh
00056       p = h* (s-k*qh) - k*rh
00057    30 IF (.NOT. (p.LT.0)) GO TO 40
00058       p = p + m
00059       GO TO 30
00060 
00061    40 GO TO 60
00062 
00063    50 p = 0
00064 C
00065 C     P = (A2*S*H)MOD M
00066 C
00067    60 IF (.NOT. (a1.NE.0)) GO TO 90
00068       q = m/a1
00069       k = s/q
00070       p = p - k* (m-a1*q)
00071       IF (p.GT.0) p = p - m
00072       p = p + a1* (s-k*q)
00073    70 IF (.NOT. (p.LT.0)) GO TO 80
00074       p = p + m
00075       GO TO 70
00076 
00077    80 CONTINUE
00078    90 k = p/qh
00079 C
00080 C     P = ((A2*H + A1)*S)MOD M
00081 C
00082       p = h* (p-k*qh) - k*rh
00083   100 IF (.NOT. (p.LT.0)) GO TO 110
00084       p = p + m
00085       GO TO 100
00086 
00087   110 CONTINUE
00088   120 IF (.NOT. (a0.NE.0)) GO TO 150
00089 C
00090 C     P = ((A2*H + A1)*H*S)MOD M
00091 C
00092       q = m/a0
00093       k = s/q
00094       p = p - k* (m-a0*q)
00095       IF (p.GT.0) p = p - m
00096       p = p + a0* (s-k*q)
00097   130 IF (.NOT. (p.LT.0)) GO TO 140
00098       p = p + m
00099       GO TO 130
00100 
00101   140 CONTINUE
00102   150 mltmod = p
00103 C
00104       RETURN
00105 
00106       END