BeamHaloGenerator/src/ranlux.f

0001       SUBROUTINE RANLUX(RVEC,LENV)
0002 C         Subtract-and-borrow random number generator proposed by
0003 C         Marsaglia and Zaman, implemented by F. James with the name
0004 C         RCARRY in 1991, and later improved by Martin Luescher
0005 C         in 1993 to produce "Luxury Pseudorandom Numbers".
0006 C     Fortran 77 coded by F. James, 1993
0007 C
0008 C   LUXURY LEVELS.
0009 C   ------ ------      The available luxury levels are:
0010 C
0011 C  level 0  (p=24): equivalent to the original RCARRY of Marsaglia
0012 C           and Zaman, very long period, but fails many tests.
0013 C  level 1  (p=48): considerable improvement in quality over level 0,
0014 C           now passes the gap test, but still fails spectral test.
0015 C  level 2  (p=97): passes all known tests, but theoretically still
0016 C           defective.
0017 C  level 3  (p=223): DEFAULT VALUE.  Any theoretically possible
0018 C           correlations have very small chance of being observed.
0019 C  level 4  (p=389): highest possible luxury, all 24 bits chaotic.
0020 C
0021 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0022 C!!!  Calling sequences for RANLUX:                                  ++
0023 C!!!      CALL RANLUX (RVEC, LEN)   returns a vector RVEC of LEN     ++
0024 C!!!                   32-bit random floating point numbers between  ++
0025 C!!!                   zero (not included) and one (also not incl.). ++
0026 C!!!      CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from  ++
0027 C!!!               one 32-bit integer INT and sets Luxury Level LUX  ++
0028 C!!!               which is integer between zero and MAXLEV, or if   ++
0029 C!!!               LUX .GT. 24, it sets p=LUX directly.  K1 and K2   ++
0030 C!!!               should be set to zero unless restarting at a break++
0031 C!!!               point given by output of RLUXAT (see RLUXAT).     ++
0032 C!!!      CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++
0033 C!!!               which can be used to restart the RANLUX generator ++
0034 C!!!               at the current point by calling RLUXGO.  K1 and K2++
0035 C!!!               specify how many numbers were generated since the ++
0036 C!!!               initialization with LUX and INT.  The restarting  ++
0037 C!!!               skips over  K1+K2*E9   numbers, so it can be long.++
0038 C!!!   A more efficient but less convenient way of restarting is by: ++
0039 C!!!      CALL RLUXIN(ISVEC)    restarts the generator from vector   ++
0040 C!!!                   ISVEC of 25 32-bit integers (see RLUXUT)      ++
0041 C!!!      CALL RLUXUT(ISVEC)    outputs the current values of the 25 ++
0042 C!!!                 32-bit integer seeds, to be used for restarting ++
0043 C!!!      ISVEC must be dimensioned 25 in the calling program        ++
0044 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0045       DIMENSION RVEC(LENV)
0046       DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25)
0047       PARAMETER (MAXLEV=4, LXDFLT=3)
0048       DIMENSION NDSKIP(0:MAXLEV)
0049       DIMENSION NEXT(24)
0050       PARAMETER (TWOP12=4096., IGIGA=1000000000,JSDFLT=314159265)
0051       PARAMETER (ITWO24=2**24, ICONS=2147483563)
0052       SAVE NOTYET, I24, J24, CARRY, SEEDS, TWOM24, TWOM12, LUXLEV
0053       SAVE NSKIP, NDSKIP, IN24, NEXT, KOUNT, MKOUNT, INSEED
0054       INTEGER LUXLEV
0055       LOGICAL NOTYET
0056       DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
0057       DATA I24,J24,CARRY/24,10,0./
0058 C                               default
0059 C  Luxury Level   0     1     2   *3*    4
0060       DATA NDSKIP/0,   24,   73,  199,  365 /
0061 Corresponds to p=24    48    97   223   389
0062 C     time factor 1     2     3     6    10   on slow workstation
0063 C                 1    1.5    2     3     5   on fast mainframe
0064 C
0065 C  NOTYET is .TRUE. if no initialization has been performed yet.
0066 C              Default Initialization by Multiplicative Congruential
0067       IF (NOTYET) THEN
0068          NOTYET = .FALSE.
0069          JSEED = JSDFLT
0070          INSEED = JSEED
0071          WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED
0072          LUXLEV = LXDFLT
0073          NSKIP = NDSKIP(LUXLEV)
0074          LP = NSKIP + 24
0075          IN24 = 0
0076          KOUNT = 0
0077          MKOUNT = 0
0078          WRITE(6,'(A,I2,A,I4)')  ' RANLUX DEFAULT LUXURY LEVEL =  ',
0079      +        LUXLEV,'      p =',LP
0080             TWOM24 = 1.
0081          DO 25 I= 1, 24
0082             TWOM24 = TWOM24 * 0.5
0083          K = JSEED/53668
0084          JSEED = 40014*(JSEED-K*53668) -K*12211
0085          IF (JSEED .LT. 0)  JSEED = JSEED+ICONS
0086          ISEEDS(I) = MOD(JSEED,ITWO24)
0087    25    CONTINUE
0088          TWOM12 = TWOM24 * 4096.
0089          DO 50 I= 1,24
0090          SEEDS(I) = REAL(ISEEDS(I))*TWOM24
0091          NEXT(I) = I-1
0092    50    CONTINUE
0093          NEXT(1) = 24
0094          I24 = 24
0095          J24 = 10
0096          CARRY = 0.
0097          IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
0098       ENDIF
0099 C
0100 C          The Generator proper: "Subtract-with-borrow",
0101 C          as proposed by Marsaglia and Zaman,
0102 C          Florida State University, March, 1989
0103 C
0104       DO 100 IVEC= 1, LENV
0105 *     UNI = SEEDS(J24) - SEEDS(I24) - CARRY
0106 *     IF (UNI .LT. 0.)  THEN
0107 *        UNI = UNI + 1.0
0108 *        CARRY = TWOM24
0109 *     ELSE
0110 *        CARRY = 0.
0111 *     ENDIF
0112 *     SEEDS(I24) = UNI
0113 *     I24 = NEXT(I24)
0114 *     J24 = NEXT(J24)
0115 *     RVEC(IVEC) = UNI
0116       uni = BHGPYR(idummy)
0117       rvec(ivec) = uni
0118 C  small numbers (with less than 12 "significant" bits) are "padded".
0119 *     IF (UNI .LT. TWOM12)  THEN
0120 *        RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
0121 C        and zero is forbidden in case someone takes a logarithm
0122 *        IF (RVEC(IVEC) .EQ. 0.)  RVEC(IVEC) = TWOM24*TWOM24
0123 *     ENDIF
0124 C        Skipping to luxury.  As proposed by Martin Luscher.
0125 *     IN24 = IN24 + 1
0126 *     IF (IN24 .EQ. 24)  THEN
0127 *        IN24 = 0
0128 *        KOUNT = KOUNT + NSKIP
0129 *        DO 90 ISK= 1, NSKIP
0130 *        UNI = SEEDS(J24) - SEEDS(I24) - CARRY
0131 *        IF (UNI .LT. 0.)  THEN
0132 *           UNI = UNI + 1.0
0133 *           CARRY = TWOM24
0134 *        ELSE
0135 *           CARRY = 0.
0136 *        ENDIF
0137 *        SEEDS(I24) = UNI
0138 *        I24 = NEXT(I24)
0139 *        J24 = NEXT(J24)
0140 *  90    CONTINUE
0141 *     ENDIF
0142   100 CONTINUE
0143       KOUNT = KOUNT + LENV
0144       IF (KOUNT .GE. IGIGA)  THEN
0145          MKOUNT = MKOUNT + 1
0146          KOUNT = KOUNT - IGIGA
0147       ENDIF
0148       RETURN
0149 C
0150 C           Entry to input and float integer seeds from previous run
0151       ENTRY RLUXIN(ISDEXT)
0152          NOTYET = .FALSE.
0153          TWOM24 = 1.
0154          DO 195 I= 1, 24
0155          NEXT(I) = I-1
0156   195    TWOM24 = TWOM24 * 0.5
0157          NEXT(1) = 24
0158          TWOM12 = TWOM24 * 4096.
0159       WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
0160       WRITE(6,'(5X,5I12)') ISDEXT
0161       DO 200 I= 1, 24
0162       SEEDS(I) = REAL(ISDEXT(I))*TWOM24
0163   200 CONTINUE
0164       CARRY = 0.
0165       IF (ISDEXT(25) .LT. 0)  CARRY = TWOM24
0166       ISD = IABS(ISDEXT(25))
0167       I24 = MOD(ISD,100)
0168       ISD = ISD/100
0169       J24 = MOD(ISD,100)
0170       ISD = ISD/100
0171       IN24 = MOD(ISD,100)
0172       ISD = ISD/100
0173       LUXLEV = ISD
0174         IF (LUXLEV .LE. MAXLEV) THEN
0175           NSKIP = NDSKIP(LUXLEV)
0176           WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',
0177      +                         LUXLEV
0178         ELSE  IF (LUXLEV .GE. 24) THEN
0179           NSKIP = LUXLEV - 24
0180           WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV
0181         ELSE
0182           NSKIP = NDSKIP(MAXLEV)
0183           WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV
0184           LUXLEV = MAXLEV
0185         ENDIF
0186       INSEED = -1
0187       RETURN
0188 C
0189 C                    Entry to ouput seeds as integers
0190       ENTRY RLUXUT(ISDEXT)
0191       DO 300 I= 1, 24
0192          ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12)
0193   300 CONTINUE
0194       ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV
0195       IF (CARRY .GT. 0.)  ISDEXT(25) = -ISDEXT(25)
0196       RETURN
0197 C
0198 C                    Entry to output the "convenient" restart point
0199       ENTRY RLUXAT(LOUT,INOUT,K1,K2)
0200       LOUT = LUXLEV
0201       INOUT = INSEED
0202       K1 = KOUNT
0203       K2 = MKOUNT
0204       RETURN
0205 C
0206 C                    Entry to initialize from one or three integers
0207       ENTRY RLUXGO(LUX,INS,K1,K2)
0208          IF (LUX .LT. 0) THEN
0209             LUXLEV = LXDFLT
0210          ELSE IF (LUX .LE. MAXLEV) THEN
0211             LUXLEV = LUX
0212          ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN
0213             LUXLEV = MAXLEV
0214             WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX
0215          ELSE
0216             LUXLEV = LUX
0217             DO 310 ILX= 0, MAXLEV
0218               IF (LUX .EQ. NDSKIP(ILX)+24)  LUXLEV = ILX
0219   310       CONTINUE
0220          ENDIF
0221       IF (LUXLEV .LE. MAXLEV)  THEN
0222          NSKIP = NDSKIP(LUXLEV)
0223          WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :',
0224      +        LUXLEV,'     P=', NSKIP+24
0225       ELSE
0226           NSKIP = LUXLEV - 24
0227           WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV
0228       ENDIF
0229       IN24 = 0
0230       IF (INS .LT. 0)  WRITE (6,'(A)')
0231      +   ' Illegal initialization by RLUXGO, negative input seed'
0232       IF (INS .GT. 0)  THEN
0233         JSEED = INS
0234         WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',
0235      +      JSEED, K1,K2
0236       ELSE
0237         JSEED = JSDFLT
0238         WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
0239       ENDIF
0240       INSEED = JSEED
0241       NOTYET = .FALSE.
0242       TWOM24 = 1.
0243          DO 325 I= 1, 24
0244            TWOM24 = TWOM24 * 0.5
0245          K = JSEED/53668
0246          JSEED = 40014*(JSEED-K*53668) -K*12211
0247          IF (JSEED .LT. 0)  JSEED = JSEED+ICONS
0248          ISEEDS(I) = MOD(JSEED,ITWO24)
0249   325    CONTINUE
0250       TWOM12 = TWOM24 * 4096.
0251          DO 350 I= 1,24
0252          SEEDS(I) = REAL(ISEEDS(I))*TWOM24
0253          NEXT(I) = I-1
0254   350    CONTINUE
0255       NEXT(1) = 24
0256       I24 = 24
0257       J24 = 10
0258       CARRY = 0.
0259       IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
0260 C        If restarting at a break point, skip K1 + IGIGA*K2
0261 C        Note that this is the number of numbers delivered to
0262 C        the user PLUS the number skipped (if luxury .GT. 0).
0263       KOUNT = K1
0264       MKOUNT = K2
0265       IF (K1+K2 .NE. 0)  THEN
0266         DO 500 IOUTER= 1, K2+1
0267           INNER = IGIGA
0268           IF (IOUTER .EQ. K2+1)  INNER = K1
0269           DO 450 ISK= 1, INNER
0270             UNI = SEEDS(J24) - SEEDS(I24) - CARRY
0271             IF (UNI .LT. 0.)  THEN
0272                UNI = UNI + 1.0
0273                CARRY = TWOM24
0274             ELSE
0275                CARRY = 0.
0276             ENDIF
0277             SEEDS(I24) = UNI
0278             I24 = NEXT(I24)
0279             J24 = NEXT(J24)
0280   450     CONTINUE
0281   500   CONTINUE
0282 C         Get the right value of IN24 by direct calculation
0283         IN24 = MOD(KOUNT, NSKIP+24)
0284         IF (MKOUNT .GT. 0)  THEN
0285            IZIP = MOD(IGIGA, NSKIP+24)
0286            IZIP2 = MKOUNT*IZIP + IN24
0287            IN24 = MOD(IZIP2, NSKIP+24)
0288         ENDIF
0289 C       Now IN24 had better be between zero and 23 inclusive
0290         IF (IN24 .GT. 23) THEN
0291            WRITE (6,'(A/A,3I11,A,I5)')
0292      +    '  Error in RESTARTING with RLUXGO:','  The values', INS,
0293      +     K1, K2, ' cannot occur at luxury level', LUXLEV
0294            IN24 = 0
0295         ENDIF
0296       ENDIF
0297       RETURN
0298       END