[u/mrichter/AliRoot.git] / MEVSIM / ranlux.F

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