[u/mrichter/AliRoot.git] / HBTP / ranlux.f

C*
C* $Id$
C*
C#include "gen/pilot.h"
      SUBROUTINE RANLUX(RVEC,LENV)
C         Subtract-and-borrow random number generator proposed by
C         Marsaglia and Zaman, implemented by F. James with the name
C         RCARRY in 1991, and later improved by Martin Luescher
C         in 1993 to produce "Luxury Pseudorandom Numbers".
C     Fortran 77 coded by F. James, 1993
C
C   LUXURY LEVELS.
C   ------ ------      The available luxury levels are:
C
C  level 0  (p=24): equivalent to the original RCARRY of Marsaglia
C           and Zaman, very long period, but fails many tests.
C  level 1  (p=48): considerable improvement in quality over level 0,
C           now passes the gap test, but still fails spectral test.
C  level 2  (p=97): passes all known tests, but theoretically still
C           defective.
C  level 3  (p=223): DEFAULT VALUE.  Any theoretically possible
C           correlations have very small chance of being observed.
C  level 4  (p=389): highest possible luxury, all 24 bits chaotic.
C
C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C!!!  Calling sequences for RANLUX:                                  ++
C!!!      CALL RANLUX (RVEC, LEN)   returns a vector RVEC of LEN     ++
C!!!                   32-bit random floating point numbers between  ++
C!!!                   zero (not included) and one (also not incl.). ++
C!!!      CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from  ++
C!!!               one 32-bit integer INT and sets Luxury Level LUX  ++
C!!!               which is integer between zero and MAXLEV, or if   ++
C!!!               LUX .GT. 24, it sets p=LUX directly.  K1 and K2   ++
C!!!               should be set to zero unless restarting at a break++ 
C!!!               point given by output of RLUXAT (see RLUXAT).     ++
C!!!      CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++
C!!!               which can be used to restart the RANLUX generator ++
C!!!               at the current point by calling RLUXGO.  K1 and K2++
C!!!               specify how many numbers were generated since the ++
C!!!               initialization with LUX and INT.  The restarting  ++
C!!!               skips over  K1+K2*E9   numbers, so it can be long.++
C!!!   A more efficient but less convenient way of restarting is by: ++
C!!!      CALL RLUXIN(ISVEC)    restarts the generator from vector   ++
C!!!                   ISVEC of 25 32-bit integers (see RLUXUT)      ++
C!!!      CALL RLUXUT(ISVEC)    outputs the current values of the 25 ++
C!!!                 32-bit integer seeds, to be used for restarting ++
C!!!      ISVEC must be dimensioned 25 in the calling program        ++
C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      DIMENSION RVEC(LENV)
      DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25)
      PARAMETER (MAXLEV=4, LXDFLT=3)
      DIMENSION NDSKIP(0:MAXLEV)
      DIMENSION NEXT(24)
      PARAMETER (TWOP12=4096., IGIGA=1000000000,JSDFLT=314159265)
      PARAMETER (ITWO24=2**24, ICONS=2147483563)
      SAVE NOTYET, I24, J24, CARRY, SEEDS, TWOM24, TWOM12, LUXLEV
      SAVE NSKIP, NDSKIP, IN24, NEXT, KOUNT, MKOUNT, INSEED
      INTEGER LUXLEV
      LOGICAL NOTYET
      DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
      DATA I24,J24,CARRY/24,10,0./
C                               default
C  Luxury Level   0     1     2   *3*    4
      DATA NDSKIP/0,   24,   73,  199,  365 /
Corresponds to p=24    48    97   223   389
C     time factor 1     2     3     6    10   on slow workstation
C                 1    1.5    2     3     5   on fast mainframe
C
C  NOTYET is .TRUE. if no initialization has been performed yet.
C              Default Initialization by Multiplicative Congruential
      IF (NOTYET) THEN
         NOTYET = .FALSE.
         JSEED = JSDFLT  
         INSEED = JSEED
         WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED
         LUXLEV = LXDFLT
         NSKIP = NDSKIP(LUXLEV)
         LP = NSKIP + 24
         IN24 = 0
         KOUNT = 0
         MKOUNT = 0
         WRITE(6,'(A,I2,A,I4)')  ' RANLUX DEFAULT LUXURY LEVEL =  ',
     +        LUXLEV,'      p =',LP
            TWOM24 = 1.
         DO 25 I= 1, 24
            TWOM24 = TWOM24 * 0.5
         K = JSEED/53668
         JSEED = 40014*(JSEED-K*53668) -K*12211
         IF (JSEED .LT. 0)  JSEED = JSEED+ICONS
         ISEEDS(I) = MOD(JSEED,ITWO24)
   25    CONTINUE
         TWOM12 = TWOM24 * 4096.
         DO 50 I= 1,24
         SEEDS(I) = REAL(ISEEDS(I))*TWOM24
         NEXT(I) = I-1
   50    CONTINUE
         NEXT(1) = 24
         I24 = 24
         J24 = 10
         CARRY = 0.
         IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
      ENDIF
C
C          The Generator proper: "Subtract-with-borrow",
C          as proposed by Marsaglia and Zaman,
C          Florida State University, March, 1989
C
      DO 100 IVEC= 1, LENV
      UNI = SEEDS(J24) - SEEDS(I24) - CARRY 
      IF (UNI .LT. 0.)  THEN
         UNI = UNI + 1.0
         CARRY = TWOM24
      ELSE
         CARRY = 0.
      ENDIF
      SEEDS(I24) = UNI
      I24 = NEXT(I24)
      J24 = NEXT(J24)
      RVEC(IVEC) = UNI
C  small numbers (with less than 12 "significant" bits) are "padded".
      IF (UNI .LT. TWOM12)  THEN
         RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
C        and zero is forbidden in case someone takes a logarithm
         IF (RVEC(IVEC) .EQ. 0.)  RVEC(IVEC) = TWOM24*TWOM24
      ENDIF
C        Skipping to luxury.  As proposed by Martin Luscher.
      IN24 = IN24 + 1
      IF (IN24 .EQ. 24)  THEN
         IN24 = 0
         KOUNT = KOUNT + NSKIP
         DO 90 ISK= 1, NSKIP
         UNI = SEEDS(J24) - SEEDS(I24) - CARRY
         IF (UNI .LT. 0.)  THEN
            UNI = UNI + 1.0
            CARRY = TWOM24
         ELSE
            CARRY = 0.
         ENDIF
         SEEDS(I24) = UNI
         I24 = NEXT(I24)
         J24 = NEXT(J24)
   90    CONTINUE
      ENDIF
  100 CONTINUE
      KOUNT = KOUNT + LENV
      IF (KOUNT .GE. IGIGA)  THEN
         MKOUNT = MKOUNT + 1
         KOUNT = KOUNT - IGIGA
      ENDIF
      RETURN
C
C           Entry to input and float integer seeds from previous run
      ENTRY RLUXIN(ISDEXT)
         NOTYET = .FALSE.
         TWOM24 = 1.
         DO 195 I= 1, 24
         NEXT(I) = I-1
  195    TWOM24 = TWOM24 * 0.5
         NEXT(1) = 24
         TWOM12 = TWOM24 * 4096.
      WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
      WRITE(6,'(5X,5I12)') ISDEXT
      DO 200 I= 1, 24
      SEEDS(I) = REAL(ISDEXT(I))*TWOM24
  200 CONTINUE
      CARRY = 0.
      IF (ISDEXT(25) .LT. 0)  CARRY = TWOM24
      ISD = IABS(ISDEXT(25))
      I24 = MOD(ISD,100)
      ISD = ISD/100
      J24 = MOD(ISD,100)
      ISD = ISD/100
      IN24 = MOD(ISD,100)
      ISD = ISD/100
      LUXLEV = ISD
        IF (LUXLEV .LE. MAXLEV) THEN
          NSKIP = NDSKIP(LUXLEV)
          WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',
     +                         LUXLEV
        ELSE  IF (LUXLEV .GE. 24) THEN
          NSKIP = LUXLEV - 24
          WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV
        ELSE
          NSKIP = NDSKIP(MAXLEV)
          WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV
          LUXLEV = MAXLEV
        ENDIF
      INSEED = -1
      RETURN
C
C                    Entry to ouput seeds as integers
      ENTRY RLUXUT(ISDEXT)
      DO 300 I= 1, 24
         ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12)
  300 CONTINUE
      ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV
      IF (CARRY .GT. 0.)  ISDEXT(25) = -ISDEXT(25)
      RETURN
C
C                    Entry to output the "convenient" restart point
      ENTRY RLUXAT(LOUT,INOUT,K1,K2)
      LOUT = LUXLEV
      INOUT = INSEED
      K1 = KOUNT
      K2 = MKOUNT
      RETURN
C
C                    Entry to initialize from one or three integers
      ENTRY RLUXGO(LUX,INS,K1,K2)
         IF (LUX .LT. 0) THEN
            LUXLEV = LXDFLT
         ELSE IF (LUX .LE. MAXLEV) THEN
            LUXLEV = LUX
         ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN
            LUXLEV = MAXLEV
            WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX
         ELSE
            LUXLEV = LUX
            DO 310 ILX= 0, MAXLEV
              IF (LUX .EQ. NDSKIP(ILX)+24)  LUXLEV = ILX
  310       CONTINUE
         ENDIF
      IF (LUXLEV .LE. MAXLEV)  THEN
         NSKIP = NDSKIP(LUXLEV)
         WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :',
     +        LUXLEV,'     P=', NSKIP+24
      ELSE
          NSKIP = LUXLEV - 24
          WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV
      ENDIF
      IN24 = 0
      IF (INS .LT. 0)  WRITE (6,'(A)')   
     +   ' Illegal initialization by RLUXGO, negative input seed'
      IF (INS .GT. 0)  THEN
        JSEED = INS
        WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',
     +      JSEED, K1,K2
      ELSE
        JSEED = JSDFLT
        WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
      ENDIF
      INSEED = JSEED
      NOTYET = .FALSE.
      TWOM24 = 1.
         DO 325 I= 1, 24
           TWOM24 = TWOM24 * 0.5
         K = JSEED/53668
         JSEED = 40014*(JSEED-K*53668) -K*12211
         IF (JSEED .LT. 0)  JSEED = JSEED+ICONS
         ISEEDS(I) = MOD(JSEED,ITWO24)
  325    CONTINUE
      TWOM12 = TWOM24 * 4096.
         DO 350 I= 1,24
         SEEDS(I) = REAL(ISEEDS(I))*TWOM24
         NEXT(I) = I-1
  350    CONTINUE
      NEXT(1) = 24
      I24 = 24
      J24 = 10
      CARRY = 0.
      IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
C        If restarting at a break point, skip K1 + IGIGA*K2
C        Note that this is the number of numbers delivered to
C        the user PLUS the number skipped (if luxury .GT. 0).
      KOUNT = K1
      MKOUNT = K2
      IF (K1+K2 .NE. 0)  THEN
        DO 500 IOUTER= 1, K2+1
          INNER = IGIGA
          IF (IOUTER .EQ. K2+1)  INNER = K1
          DO 450 ISK= 1, INNER
            UNI = SEEDS(J24) - SEEDS(I24) - CARRY 
            IF (UNI .LT. 0.)  THEN
               UNI = UNI + 1.0
               CARRY = TWOM24
            ELSE
               CARRY = 0.
            ENDIF
            SEEDS(I24) = UNI
            I24 = NEXT(I24)
            J24 = NEXT(J24)
  450     CONTINUE
  500   CONTINUE
C         Get the right value of IN24 by direct calculation
        IN24 = MOD(KOUNT, NSKIP+24)
        IF (MKOUNT .GT. 0)  THEN
           IZIP = MOD(IGIGA, NSKIP+24)
           IZIP2 = MKOUNT*IZIP + IN24
           IN24 = MOD(IZIP2, NSKIP+24)
        ENDIF
C       Now IN24 had better be between zero and 23 inclusive
        IF (IN24 .GT. 23) THEN
           WRITE (6,'(A/A,3I11,A,I5)')  
     +    '  Error in RESTARTING with RLUXGO:','  The values', INS,
     +     K1, K2, ' cannot occur at luxury level', LUXLEV
           IN24 = 0
        ENDIF
      ENDIF
      RETURN
      END
Commit	Line	Data
18448239	1	C*
	2	C* $Id$
	3	C*
18448239	4	C#include "gen/pilot.h"
	5	SUBROUTINE RANLUX(RVEC,LENV)
	6	C Subtract-and-borrow random number generator proposed by
	7	C Marsaglia and Zaman, implemented by F. James with the name
	8	C RCARRY in 1991, and later improved by Martin Luescher
	9	C in 1993 to produce "Luxury Pseudorandom Numbers".
	10	C Fortran 77 coded by F. James, 1993
	11	C
	12	C LUXURY LEVELS.
	13	C ------ ------ The available luxury levels are:
	14	C
	15	C level 0 (p=24): equivalent to the original RCARRY of Marsaglia
	16	C and Zaman, very long period, but fails many tests.
	17	C level 1 (p=48): considerable improvement in quality over level 0,
	18	C now passes the gap test, but still fails spectral test.
	19	C level 2 (p=97): passes all known tests, but theoretically still
	20	C defective.
	21	C level 3 (p=223): DEFAULT VALUE. Any theoretically possible
	22	C correlations have very small chance of being observed.
	23	C level 4 (p=389): highest possible luxury, all 24 bits chaotic.
	24	C
	25	C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	26	C!!! Calling sequences for RANLUX: ++
	27	C!!! CALL RANLUX (RVEC, LEN) returns a vector RVEC of LEN ++
	28	C!!! 32-bit random floating point numbers between ++
	29	C!!! zero (not included) and one (also not incl.). ++
	30	C!!! CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from ++
	31	C!!! one 32-bit integer INT and sets Luxury Level LUX ++
	32	C!!! which is integer between zero and MAXLEV, or if ++
	33	C!!! LUX .GT. 24, it sets p=LUX directly. K1 and K2 ++
	34	C!!! should be set to zero unless restarting at a break++
	35	C!!! point given by output of RLUXAT (see RLUXAT). ++
	36	C!!! CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++
	37	C!!! which can be used to restart the RANLUX generator ++
	38	C!!! at the current point by calling RLUXGO. K1 and K2++
	39	C!!! specify how many numbers were generated since the ++
	40	C!!! initialization with LUX and INT. The restarting ++
	41	C!!! skips over K1+K2*E9 numbers, so it can be long.++
	42	C!!! A more efficient but less convenient way of restarting is by: ++
	43	C!!! CALL RLUXIN(ISVEC) restarts the generator from vector ++
	44	C!!! ISVEC of 25 32-bit integers (see RLUXUT) ++
	45	C!!! CALL RLUXUT(ISVEC) outputs the current values of the 25 ++
	46	C!!! 32-bit integer seeds, to be used for restarting ++
	47	C!!! ISVEC must be dimensioned 25 in the calling program ++
	48	C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	49	DIMENSION RVEC(LENV)
	50	DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25)
	51	PARAMETER (MAXLEV=4, LXDFLT=3)
	52	DIMENSION NDSKIP(0:MAXLEV)
	53	DIMENSION NEXT(24)
	54	PARAMETER (TWOP12=4096., IGIGA=1000000000,JSDFLT=314159265)
	55	PARAMETER (ITWO24=2**24, ICONS=2147483563)
	56	SAVE NOTYET, I24, J24, CARRY, SEEDS, TWOM24, TWOM12, LUXLEV
	57	SAVE NSKIP, NDSKIP, IN24, NEXT, KOUNT, MKOUNT, INSEED
	58	INTEGER LUXLEV
	59	LOGICAL NOTYET
	60	DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
	61	DATA I24,J24,CARRY/24,10,0./
	62	C default
	63	C Luxury Level 0 1 2 3 4
	64	DATA NDSKIP/0, 24, 73, 199, 365 /
	65	Corresponds to p=24 48 97 223 389
	66	C time factor 1 2 3 6 10 on slow workstation
	67	C 1 1.5 2 3 5 on fast mainframe
68	C
69	C NOTYET is .TRUE. if no initialization has been performed yet.
70	C Default Initialization by Multiplicative Congruential
71	IF (NOTYET) THEN
72	NOTYET = .FALSE.
73	JSEED = JSDFLT
74	INSEED = JSEED
75	WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED
76	LUXLEV = LXDFLT
77	NSKIP = NDSKIP(LUXLEV)
78	LP = NSKIP + 24
79	IN24 = 0
80	KOUNT = 0
81	MKOUNT = 0
82	WRITE(6,'(A,I2,A,I4)') ' RANLUX DEFAULT LUXURY LEVEL = ',
83	+ LUXLEV,' p =',LP
84	TWOM24 = 1.
85	DO 25 I= 1, 24
86	TWOM24 = TWOM24 * 0.5
87	K = JSEED/53668
88	JSEED = 40014(JSEED-K53668) -K*12211
89	IF (JSEED .LT. 0) JSEED = JSEED+ICONS
90	ISEEDS(I) = MOD(JSEED,ITWO24)
91	25 CONTINUE
92	TWOM12 = TWOM24 * 4096.
93	DO 50 I= 1,24
94	SEEDS(I) = REAL(ISEEDS(I))*TWOM24
95	NEXT(I) = I-1
96	50 CONTINUE
97	NEXT(1) = 24
98	I24 = 24
99	J24 = 10
100	CARRY = 0.
101	IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
102	ENDIF
103	C
104	C The Generator proper: "Subtract-with-borrow",
105	C as proposed by Marsaglia and Zaman,
106	C Florida State University, March, 1989
107	C
108	DO 100 IVEC= 1, LENV
109	UNI = SEEDS(J24) - SEEDS(I24) - CARRY
110	IF (UNI .LT. 0.) THEN
111	UNI = UNI + 1.0
112	CARRY = TWOM24
113	ELSE
114	CARRY = 0.
115	ENDIF
116	SEEDS(I24) = UNI
117	I24 = NEXT(I24)
118	J24 = NEXT(J24)
119	RVEC(IVEC) = UNI
120	C small numbers (with less than 12 "significant" bits) are "padded".
121	IF (UNI .LT. TWOM12) THEN
122	RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
123	C and zero is forbidden in case someone takes a logarithm
124	IF (RVEC(IVEC) .EQ. 0.) RVEC(IVEC) = TWOM24*TWOM24
125	ENDIF
126	C Skipping to luxury. As proposed by Martin Luscher.
127	IN24 = IN24 + 1
128	IF (IN24 .EQ. 24) THEN
129	IN24 = 0
130	KOUNT = KOUNT + NSKIP
131	DO 90 ISK= 1, NSKIP
132	UNI = SEEDS(J24) - SEEDS(I24) - CARRY
133	IF (UNI .LT. 0.) THEN
134	UNI = UNI + 1.0
135	CARRY = TWOM24
136	ELSE
137	CARRY = 0.
138	ENDIF
139	SEEDS(I24) = UNI
140	I24 = NEXT(I24)
141	J24 = NEXT(J24)
142	90 CONTINUE
143	ENDIF
144	100 CONTINUE
145	KOUNT = KOUNT + LENV
146	IF (KOUNT .GE. IGIGA) THEN
147	MKOUNT = MKOUNT + 1
148	KOUNT = KOUNT - IGIGA
149	ENDIF
150	RETURN
151	C
152	C Entry to input and float integer seeds from previous run
153	ENTRY RLUXIN(ISDEXT)
154	NOTYET = .FALSE.
155	TWOM24 = 1.
156	DO 195 I= 1, 24
157	NEXT(I) = I-1
158	195 TWOM24 = TWOM24 * 0.5
159	NEXT(1) = 24
160	TWOM12 = TWOM24 * 4096.
161	WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
162	WRITE(6,'(5X,5I12)') ISDEXT
163	DO 200 I= 1, 24
164	SEEDS(I) = REAL(ISDEXT(I))*TWOM24
165	200 CONTINUE
166	CARRY = 0.
167	IF (ISDEXT(25) .LT. 0) CARRY = TWOM24
168	ISD = IABS(ISDEXT(25))
169	I24 = MOD(ISD,100)
170	ISD = ISD/100
171	J24 = MOD(ISD,100)
172	ISD = ISD/100
173	IN24 = MOD(ISD,100)
174	ISD = ISD/100
175	LUXLEV = ISD
176	IF (LUXLEV .LE. MAXLEV) THEN
177	NSKIP = NDSKIP(LUXLEV)
178	WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',
179	+ LUXLEV
180	ELSE IF (LUXLEV .GE. 24) THEN
181	NSKIP = LUXLEV - 24
182	WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV
183	ELSE
184	NSKIP = NDSKIP(MAXLEV)
185	WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV
186	LUXLEV = MAXLEV
187	ENDIF
188	INSEED = -1
189	RETURN
190	C
191	C Entry to ouput seeds as integers
192	ENTRY RLUXUT(ISDEXT)
193	DO 300 I= 1, 24
194	ISDEXT(I) = INT(SEEDS(I)TWOP12TWOP12)
195	300 CONTINUE
196	ISDEXT(25) = I24 + 100J24 + 10000IN24 + 1000000*LUXLEV
197	IF (CARRY .GT. 0.) ISDEXT(25) = -ISDEXT(25)
198	RETURN
199	C
200	C Entry to output the "convenient" restart point
201	ENTRY RLUXAT(LOUT,INOUT,K1,K2)
202	LOUT = LUXLEV
203	INOUT = INSEED
204	K1 = KOUNT
205	K2 = MKOUNT
206	RETURN
207	C
208	C Entry to initialize from one or three integers
209	ENTRY RLUXGO(LUX,INS,K1,K2)
210	IF (LUX .LT. 0) THEN
211	LUXLEV = LXDFLT
212	ELSE IF (LUX .LE. MAXLEV) THEN
213	LUXLEV = LUX
214	ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN
215	LUXLEV = MAXLEV
216	WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX
217	ELSE
218	LUXLEV = LUX
219	DO 310 ILX= 0, MAXLEV
220	IF (LUX .EQ. NDSKIP(ILX)+24) LUXLEV = ILX
221	310 CONTINUE
222	ENDIF
223	IF (LUXLEV .LE. MAXLEV) THEN
224	NSKIP = NDSKIP(LUXLEV)
225	WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :',
226	+ LUXLEV,' P=', NSKIP+24
227	ELSE
228	NSKIP = LUXLEV - 24
229	WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV
230	ENDIF
231	IN24 = 0
232	IF (INS .LT. 0) WRITE (6,'(A)')
233	+ ' Illegal initialization by RLUXGO, negative input seed'
234	IF (INS .GT. 0) THEN
235	JSEED = INS
236	WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',
237	+ JSEED, K1,K2
238	ELSE
239	JSEED = JSDFLT
240	WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
241	ENDIF
242	INSEED = JSEED
243	NOTYET = .FALSE.
244	TWOM24 = 1.
245	DO 325 I= 1, 24
246	TWOM24 = TWOM24 * 0.5
247	K = JSEED/53668
248	JSEED = 40014(JSEED-K53668) -K*12211
249	IF (JSEED .LT. 0) JSEED = JSEED+ICONS
250	ISEEDS(I) = MOD(JSEED,ITWO24)
251	325 CONTINUE
252	TWOM12 = TWOM24 * 4096.
253	DO 350 I= 1,24
254	SEEDS(I) = REAL(ISEEDS(I))*TWOM24
255	NEXT(I) = I-1
256	350 CONTINUE
257	NEXT(1) = 24
258	I24 = 24
259	J24 = 10
260	CARRY = 0.
261	IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
262	C If restarting at a break point, skip K1 + IGIGA*K2
263	C Note that this is the number of numbers delivered to
264	C the user PLUS the number skipped (if luxury .GT. 0).
265	KOUNT = K1
266	MKOUNT = K2
267	IF (K1+K2 .NE. 0) THEN
268	DO 500 IOUTER= 1, K2+1
269	INNER = IGIGA
270	IF (IOUTER .EQ. K2+1) INNER = K1
271	DO 450 ISK= 1, INNER
272	UNI = SEEDS(J24) - SEEDS(I24) - CARRY
273	IF (UNI .LT. 0.) THEN
274	UNI = UNI + 1.0
275	CARRY = TWOM24
276	ELSE
277	CARRY = 0.
278	ENDIF
279	SEEDS(I24) = UNI
280	I24 = NEXT(I24)
281	J24 = NEXT(J24)
282	450 CONTINUE
283	500 CONTINUE
284	C Get the right value of IN24 by direct calculation
285	IN24 = MOD(KOUNT, NSKIP+24)
286	IF (MKOUNT .GT. 0) THEN
287	IZIP = MOD(IGIGA, NSKIP+24)
288	IZIP2 = MKOUNT*IZIP + IN24
289	IN24 = MOD(IZIP2, NSKIP+24)
290	ENDIF
291	C Now IN24 had better be between zero and 23 inclusive
292	IF (IN24 .GT. 23) THEN
293	WRITE (6,'(A/A,3I11,A,I5)')
294	+ ' Error in RESTARTING with RLUXGO:',' The values', INS,
295	+ K1, K2, ' cannot occur at luxury level', LUXLEV
296	IN24 = 0
297	ENDIF
298	ENDIF
299	RETURN
300	END