MINICERN/mathlib/gen/v/rm48.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.2  1996/12/12 16:32:06  cernlib
   6 * Variables ONE and ZERO added to SAVE statement, courtesy R.Veenhof
   7 *
   8 * Revision 1.1.1.1  1996/04/01 15:02:55  mclareni
   9 * Mathlib gen
  10 *
  11 *
  12 #include "gen/pilot.h"
  13       SUBROUTINE RM48(RVEC,LENV)
  14 C     Double-precision version of
  15 C Universal random number generator proposed by Marsaglia and Zaman
  16 C in report FSU-SCRI-87-50
  17 C        based on RANMAR, modified by F. James, to generate vectors
  18 C        of pseudorandom numbers RVEC of length LENV, where the numbers
  19 C        in RVEC are numbers with at least 48-bit mantissas.
  20 C   Input and output entry points: RM48IN, RM48UT.
  21 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  22 C!!!  Calling sequences for RM48:                                    ++
  23 C!!!      CALL RM48 (RVEC, LEN)     returns a vector RVEC of LEN     ++
  24 C!!!                   64-bit random floating point numbers between  ++
  25 C!!!                   zero and one.                                 ++
  26 C!!!      CALL RM48IN(I1,N1,N2)   initializes the generator from one ++
  27 C!!!                   64-bit integer I1, and number counts N1,N2    ++
  28 C!!!                  (for initializing, set N1=N2=0, but to restart ++
  29 C!!!                    a previously generated sequence, use values  ++
  30 C!!!                    output by RM48UT)                            ++
  31 C!!!      CALL RM48UT(I1,N1,N2)   outputs the value of the original  ++
  32 C!!!                  seed and the two number counts, to be used     ++
  33 C!!!                  for restarting by initializing to I1 and       ++
  34 C!!!                  skipping N2*100000000+N1 numbers.              ++
  35 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  36 C for 32-bit machines, use IMPLICIT DOUBLE PRECISION
  37       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  38       DIMENSION RVEC(*)
  39       COMMON/R48ST1/U(97),C,I97,J97
  40       PARAMETER (MODCNS=1000000000)
  41       SAVE CD, CM, TWOM24, NTOT, NTOT2, IJKL,TWOM49, ONE, ZERO
  42       DATA NTOT,NTOT2,IJKL/-1,0,0/
  43 C
  44       IF (NTOT .GE. 0)  GO TO 50
  45 C
  46 C        Default initialization. User has called RM48 without RM48IN.
  47       IJKL = 54217137
  48       NTOT = 0
  49       NTOT2 = 0
  50       KALLED = 0
  51       GO TO 1
  52 C
  53       ENTRY      RM48IN(IJKLIN, NTOTIN,NTOT2N)
  54 C         Initializing routine for RM48, may be called before
  55 C         generating pseudorandom numbers with RM48.   The input
  56 C         values should be in the ranges:  0<=IJKLIN<=900 OOO OOO
  57 C                                          0<=NTOTIN<=999 999 999
  58 C                                          0<=NTOT2N<<999 999 999!
  59 C To get the standard values in Marsaglia's paper, IJKLIN=54217137
  60 C                                            NTOTIN,NTOT2N=0
  61       IJKL = IJKLIN
  62       NTOT = MAX(NTOTIN,0)
  63       NTOT2= MAX(NTOT2N,0)
  64       KALLED = 1
  65 C          always come here to initialize
  66     1 CONTINUE
  67       IJ = IJKL/30082
  68       KL = IJKL - 30082*IJ
  69       I = MOD(IJ/177, 177) + 2
  70       J = MOD(IJ, 177)     + 2
  71       K = MOD(KL/169, 178) + 1
  72       L = MOD(KL, 169)
  73       WRITE(6,'(A,I10,2X,2I10)') ' RM48 INITIALIZED:',IJKL,NTOT,NTOT2
  74 CCC      PRINT '(A,4I10)', '   I,J,K,L= ',I,J,K,L
  75       ONE = 1.
  76       HALF = 0.5
  77       ZERO = 0.
  78       DO 2 II= 1, 97
  79       S = 0.
  80       T = HALF
  81       DO 3 JJ= 1, 48
  82          M = MOD(MOD(I*J,179)*K, 179)
  83          I = J
  84          J = K
  85          K = M
  86          L = MOD(53*L+1, 169)
  87          IF (MOD(L*M,64) .GE. 32)  S = S+T
  88     3    T = HALF*T
  89     2 U(II) = S
  90       TWOM49 = T
  91       TWOM24 = ONE
  92       DO 4 I24= 1, 24
  93     4 TWOM24 = HALF*TWOM24
  94       C  =   362436.*TWOM24
  95       CD =  7654321.*TWOM24
  96       CM = 16777213.*TWOM24
  97       I97 = 97
  98       J97 = 33
  99 C       Complete initialization by skipping
 100 C            (NTOT2*MODCNS + NTOT) random numbers
 101       DO 45 LOOP2= 1, NTOT2+1
 102       NOW = MODCNS
 103       IF (LOOP2 .EQ. NTOT2+1)  NOW=NTOT
 104       IF (NOW .GT. 0)  THEN
 105       WRITE(6,'(A,I15)') ' RM48IN SKIPPING OVER ',NOW
 106           DO 40 IDUM = 1, NTOT
 107           UNI = U(I97)-U(J97)
 108           IF (UNI .LT. ZERO)  UNI=UNI+ONE
 109           U(I97) = UNI
 110           I97 = I97-1
 111           IF (I97 .EQ. 0)  I97=97
 112           J97 = J97-1
 113           IF (J97 .EQ. 0)  J97=97
 114           C = C - CD
 115           IF (C .LT. ZERO)  C=C+CM
 116    40     CONTINUE
 117       ENDIF
 118    45 CONTINUE
 119       IF (KALLED .EQ. 1)  RETURN
 120 C
 121 C          Normal entry to generate LENV random numbers
 122    50 CONTINUE
 123       DO 100 IVEC= 1, LENV
 124       UNI = U(I97)-U(J97)
 125       IF (UNI .LT. ZERO)  UNI=UNI+ONE
 126       U(I97) = UNI
 127       I97 = I97-1
 128       IF (I97 .EQ. 0)  I97=97
 129       J97 = J97-1
 130       IF (J97 .EQ. 0)  J97=97
 131       C = C - CD
 132       IF (C .LT. ZERO)  C=C+CM
 133       UNI = UNI-C
 134       IF (UNI .LT. ZERO) UNI=UNI+ONE
 135       RVEC(IVEC) = UNI
 136 C             Replace exact zeros by 2**-49
 137          IF (UNI .EQ. ZERO)  THEN
 138             RVEC(IVEC) = TWOM49
 139          ENDIF
 140   100 CONTINUE
 141       NTOT = NTOT + LENV
 142          IF (NTOT .GE. MODCNS)  THEN
 143          NTOT2 = NTOT2 + 1
 144          NTOT = NTOT - MODCNS
 145          ENDIF
 146       RETURN
 147 C           Entry to output current status
 148       ENTRY RM48UT(IJKLUT,NTOTUT,NTOT2T)
 149       IJKLUT = IJKL
 150       NTOTUT = NTOT
 151       NTOT2T = NTOT2
 152       RETURN
 153       END