home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / ranlib / mltmod.f < prev    next >
Text File  |  1996-09-28  |  2KB  |  107 lines
  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
  107.