Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:29:35

0001 CCC
0002 CCC from http://cernlib.web.cern.ch/cernlib/download/2005_source/src/mathlib/gen/v/rm48.F
0003 CCC
0004 *
0005 *
0006 * Revision 1.2  2008/04/09 14:11:58  marafino
0007 * Insert new and changed files to accomodate running this code as a producer
0008 * or as a source.
0009 *
0010 * Revision 1.1  2007/02/06 14:40:25  eperez
0011 * first version
0012 *
0013 * Revision 1.2  1996/12/12 16:32:06  cernlib
0014 * Variables ONE and ZERO added to SAVE statement, courtesy R.Veenhof
0015 *
0016 * Revision 1.1.1.1  1996/04/01 15:02:55  mclareni
0017 * Mathlib gen
0018 *
0019 *
0020       SUBROUTINE RM48(RVEC,LENV)
0021 C     Double-precision version of
0022 C Universal random number generator proposed by Marsaglia and Zaman
0023 C in report FSU-SCRI-87-50
0024 C        based on RANMAR, modified by F. James, to generate vectors
0025 C        of pseudorandom numbers RVEC of length LENV, where the numbers
0026 C        in RVEC are numbers with at least 48-bit mantissas.
0027 C   Input and output entry points: RM48IN, RM48UT.
0028 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0029 C!!!  Calling sequences for RM48:                                    ++
0030 C!!!      CALL RM48 (RVEC, LEN)     returns a vector RVEC of LEN     ++
0031 C!!!                   64-bit random floating point numbers between  ++
0032 C!!!                   zero and one.                                 ++
0033 C!!!      CALL RM48IN(I1,N1,N2)   initializes the generator from one ++
0034 C!!!                   64-bit integer I1, and number counts N1,N2    ++
0035 C!!!                  (for initializing, set N1=N2=0, but to restart ++
0036 C!!!                    a previously generated sequence, use values  ++ 
0037 C!!!                    output by RM48UT)                            ++ 
0038 C!!!      CALL RM48UT(I1,N1,N2)   outputs the value of the original  ++
0039 C!!!                  seed and the two number counts, to be used     ++
0040 C!!!                  for restarting by initializing to I1 and       ++  
0041 C!!!                  skipping N2*100000000+N1 numbers.              ++
0042 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0043 C for 32-bit machines, use IMPLICIT DOUBLE PRECISION
0044       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
0045       DIMENSION RVEC(*)
0046       COMMON/R48ST1/U(97),C,I97,J97
0047       PARAMETER (MODCNS=1000000000)
0048       SAVE CD, CM, TWOM24, NTOT, NTOT2, IJKL,TWOM49, ONE, ZERO
0049       DATA NTOT,NTOT2,IJKL/-1,0,0/
0050 C
0051       IF (NTOT .GE. 0)  GO TO 50
0052 C
0053 C        Default initialization. User has called RM48 without RM48IN.
0054       IJKL = 54217137
0055       NTOT = 0
0056       NTOT2 = 0
0057       KALLED = 0
0058       GO TO 1
0059 C
0060       ENTRY      RM48IN(IJKLIN, NTOTIN,NTOT2N)
0061 C         Initializing routine for RM48, may be called before
0062 C         generating pseudorandom numbers with RM48.   The input
0063 C         values should be in the ranges:  0<=IJKLIN<=900 OOO OOO
0064 C                                          0<=NTOTIN<=999 999 999
0065 C                                          0<=NTOT2N<<999 999 999!
0066 C To get the standard values in Marsaglia's paper, IJKLIN=54217137
0067 C                                            NTOTIN,NTOT2N=0
0068       IJKL = IJKLIN
0069       NTOT = MAX(NTOTIN,0)
0070       NTOT2= MAX(NTOT2N,0)
0071       KALLED = 1
0072 C          always come here to initialize
0073     1 CONTINUE
0074       IJ = IJKL/30082
0075       KL = IJKL - 30082*IJ
0076       I = MOD(IJ/177, 177) + 2
0077       J = MOD(IJ, 177)     + 2
0078       K = MOD(KL/169, 178) + 1
0079       L = MOD(KL, 169)
0080       WRITE(6,'(A,I10,2X,2I10)') ' RM48 INITIALIZED:',IJKL,NTOT,NTOT2
0081 CCC      PRINT '(A,4I10)', '   I,J,K,L= ',I,J,K,L
0082       ONE = 1.
0083       HALF = 0.5
0084       ZERO = 0.
0085       DO 2 II= 1, 97
0086       S = 0.
0087       T = HALF
0088       DO 3 JJ= 1, 48
0089          M = MOD(MOD(I*J,179)*K, 179)
0090          I = J
0091          J = K
0092          K = M
0093          L = MOD(53*L+1, 169)
0094          IF (MOD(L*M,64) .GE. 32)  S = S+T
0095     3    T = HALF*T
0096     2 U(II) = S
0097       TWOM49 = T
0098       TWOM24 = ONE
0099       DO 4 I24= 1, 24
0100     4 TWOM24 = HALF*TWOM24
0101       C  =   362436.*TWOM24
0102       CD =  7654321.*TWOM24
0103       CM = 16777213.*TWOM24
0104       I97 = 97
0105       J97 = 33
0106 C       Complete initialization by skipping
0107 C            (NTOT2*MODCNS + NTOT) random numbers
0108       DO 45 LOOP2= 1, NTOT2+1
0109       NOW = MODCNS
0110       IF (LOOP2 .EQ. NTOT2+1)  NOW=NTOT
0111       IF (NOW .GT. 0)  THEN
0112       WRITE(6,'(A,I15)') ' RM48IN SKIPPING OVER ',NOW
0113           DO 40 IDUM = 1, NTOT
0114           UNI = U(I97)-U(J97)
0115           IF (UNI .LT. ZERO)  UNI=UNI+ONE
0116           U(I97) = UNI
0117           I97 = I97-1
0118           IF (I97 .EQ. 0)  I97=97
0119           J97 = J97-1
0120           IF (J97 .EQ. 0)  J97=97
0121           C = C - CD
0122           IF (C .LT. ZERO)  C=C+CM
0123    40     CONTINUE
0124       ENDIF
0125    45 CONTINUE
0126       IF (KALLED .EQ. 1)  RETURN
0127 C
0128 C          Normal entry to generate LENV random numbers
0129 C          Modified 13 Marcch '08 to use the CLHEP engines
0130    50 CONTINUE
0131       DO 100 IVEC= 1, LENV
0132 *     UNI = U(I97)-U(J97)
0133 *     IF (UNI .LT. ZERO)  UNI=UNI+ONE
0134 *     U(I97) = UNI
0135 *     I97 = I97-1
0136 *     IF (I97 .EQ. 0)  I97=97
0137 *     J97 = J97-1
0138 *     IF (J97 .EQ. 0)  J97=97
0139 *     C = C - CD
0140 *     IF (C .LT. ZERO)  C=C+CM
0141 *     UNI = UNI-C
0142 *     IF (UNI .LT. ZERO) UNI=UNI+ONE
0143 *     RVEC(IVEC) = UNI
0144       UNI = BHGPYR(IDUMMY)
0145       RVEC(IVEC) = UNI
0146 C             Replace exact zeros by 2**-49
0147 *        IF (UNI .EQ. ZERO)  THEN
0148 *           RVEC(IVEC) = TWOM49
0149 *        ENDIF
0150   100 CONTINUE
0151       NTOT = NTOT + LENV
0152          IF (NTOT .GE. MODCNS)  THEN
0153          NTOT2 = NTOT2 + 1
0154          NTOT = NTOT - MODCNS
0155          ENDIF
0156       RETURN
0157 C           Entry to output current status
0158       ENTRY RM48UT(IJKLUT,NTOTUT,NTOT2T)
0159       IJKLUT = IJKL
0160       NTOTUT = NTOT
0161       NTOT2T = NTOT2
0162       RETURN
0163       END