[u/mrichter/AliRoot.git] / HBTP / ranlux2.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 RANLUX2(RVEC,LENV,Input_seed)
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
      Integer Input_seed,JSDFLT_set
      LOGICAL NOTYET
      DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
      DATA I24,J24,CARRY/24,10,0./
CCC  Set starting seed value:
      If(Input_seed.gt.0) then
         JSDFLT_set = Input_seed
      Else
         JSDFLT_set = JSDFLT
      End If
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_set  
         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_set
        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*
	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 RANLUX2(RVEC,LENV,Input_seed)
	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	Integer Input_seed,JSDFLT_set
69	LOGICAL NOTYET
70	DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
71	DATA I24,J24,CARRY/24,10,0./
72	CCC Set starting seed value:
73	If(Input_seed.gt.0) then
74	JSDFLT_set = Input_seed
75	Else
76	JSDFLT_set = JSDFLT
77	End If
78	C default
79	C Luxury Level 0 1 2 3 4
80	DATA NDSKIP/0, 24, 73, 199, 365 /
81	Corresponds to p=24 48 97 223 389
82	C time factor 1 2 3 6 10 on slow workstation
83	C 1 1.5 2 3 5 on fast mainframe
84	C
85	C NOTYET is .TRUE. if no initialization has been performed yet.
86	C Default Initialization by Multiplicative Congruential
87	IF (NOTYET) THEN
88	NOTYET = .FALSE.
89	JSEED = JSDFLT_set
90	INSEED = JSEED
91	WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED
92	LUXLEV = LXDFLT
93	NSKIP = NDSKIP(LUXLEV)
94	LP = NSKIP + 24
95	IN24 = 0
96	KOUNT = 0
97	MKOUNT = 0
98	WRITE(6,'(A,I2,A,I4)') ' RANLUX DEFAULT LUXURY LEVEL = ',
99	+ LUXLEV,' p =',LP
100	TWOM24 = 1.
101	DO 25 I= 1, 24
102	TWOM24 = TWOM24 * 0.5
103	K = JSEED/53668
104	JSEED = 40014(JSEED-K53668) -K*12211
105	IF (JSEED .LT. 0) JSEED = JSEED+ICONS
106	ISEEDS(I) = MOD(JSEED,ITWO24)
107	25 CONTINUE
108	TWOM12 = TWOM24 * 4096.
109	DO 50 I= 1,24
110	SEEDS(I) = REAL(ISEEDS(I))*TWOM24
111	NEXT(I) = I-1
112	50 CONTINUE
113	NEXT(1) = 24
114	I24 = 24
115	J24 = 10
116	CARRY = 0.
117	IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
118	ENDIF
119	C
120	C The Generator proper: "Subtract-with-borrow",
121	C as proposed by Marsaglia and Zaman,
122	C Florida State University, March, 1989
123	C
124	DO 100 IVEC= 1, LENV
125	UNI = SEEDS(J24) - SEEDS(I24) - CARRY
126	IF (UNI .LT. 0.) THEN
127	UNI = UNI + 1.0
128	CARRY = TWOM24
129	ELSE
130	CARRY = 0.
131	ENDIF
132	SEEDS(I24) = UNI
133	I24 = NEXT(I24)
134	J24 = NEXT(J24)
135	RVEC(IVEC) = UNI
136	C small numbers (with less than 12 "significant" bits) are "padded".
137	IF (UNI .LT. TWOM12) THEN
138	RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
139	C and zero is forbidden in case someone takes a logarithm
140	IF (RVEC(IVEC) .EQ. 0.) RVEC(IVEC) = TWOM24*TWOM24
141	ENDIF
142	C Skipping to luxury. As proposed by Martin Luscher.
143	IN24 = IN24 + 1
144	IF (IN24 .EQ. 24) THEN
145	IN24 = 0
146	KOUNT = KOUNT + NSKIP
147	DO 90 ISK= 1, NSKIP
148	UNI = SEEDS(J24) - SEEDS(I24) - CARRY
149	IF (UNI .LT. 0.) THEN
150	UNI = UNI + 1.0
151	CARRY = TWOM24
152	ELSE
153	CARRY = 0.
154	ENDIF
155	SEEDS(I24) = UNI
156	I24 = NEXT(I24)
157	J24 = NEXT(J24)
158	90 CONTINUE
159	ENDIF
160	100 CONTINUE
161	KOUNT = KOUNT + LENV
162	IF (KOUNT .GE. IGIGA) THEN
163	MKOUNT = MKOUNT + 1
164	KOUNT = KOUNT - IGIGA
165	ENDIF
166	RETURN
167	C
168	C Entry to input and float integer seeds from previous run
169	ENTRY RLUXIN(ISDEXT)
170	NOTYET = .FALSE.
171	TWOM24 = 1.
172	DO 195 I= 1, 24
173	NEXT(I) = I-1
174	195 TWOM24 = TWOM24 * 0.5
175	NEXT(1) = 24
176	TWOM12 = TWOM24 * 4096.
177	WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
178	WRITE(6,'(5X,5I12)') ISDEXT
179	DO 200 I= 1, 24
180	SEEDS(I) = REAL(ISDEXT(I))*TWOM24
181	200 CONTINUE
182	CARRY = 0.
183	IF (ISDEXT(25) .LT. 0) CARRY = TWOM24
184	ISD = IABS(ISDEXT(25))
185	I24 = MOD(ISD,100)
186	ISD = ISD/100
187	J24 = MOD(ISD,100)
188	ISD = ISD/100
189	IN24 = MOD(ISD,100)
190	ISD = ISD/100
191	LUXLEV = ISD
192	IF (LUXLEV .LE. MAXLEV) THEN
193	NSKIP = NDSKIP(LUXLEV)
194	WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',
195	+ LUXLEV
196	ELSE IF (LUXLEV .GE. 24) THEN
197	NSKIP = LUXLEV - 24
198	WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV
199	ELSE
200	NSKIP = NDSKIP(MAXLEV)
201	WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV
202	LUXLEV = MAXLEV
203	ENDIF
204	INSEED = -1
205	RETURN
206	C
207	C Entry to ouput seeds as integers
208	ENTRY RLUXUT(ISDEXT)
209	DO 300 I= 1, 24
210	ISDEXT(I) = INT(SEEDS(I)TWOP12TWOP12)
211	300 CONTINUE
212	ISDEXT(25) = I24 + 100J24 + 10000IN24 + 1000000*LUXLEV
213	IF (CARRY .GT. 0.) ISDEXT(25) = -ISDEXT(25)
214	RETURN
215	C
216	C Entry to output the "convenient" restart point
217	ENTRY RLUXAT(LOUT,INOUT,K1,K2)
218	LOUT = LUXLEV
219	INOUT = INSEED
220	K1 = KOUNT
221	K2 = MKOUNT
222	RETURN
223	C
224	C Entry to initialize from one or three integers
225	ENTRY RLUXGO(LUX,INS,K1,K2)
226	IF (LUX .LT. 0) THEN
227	LUXLEV = LXDFLT
228	ELSE IF (LUX .LE. MAXLEV) THEN
229	LUXLEV = LUX
230	ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN
231	LUXLEV = MAXLEV
232	WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX
233	ELSE
234	LUXLEV = LUX
235	DO 310 ILX= 0, MAXLEV
236	IF (LUX .EQ. NDSKIP(ILX)+24) LUXLEV = ILX
237	310 CONTINUE
238	ENDIF
239	IF (LUXLEV .LE. MAXLEV) THEN
240	NSKIP = NDSKIP(LUXLEV)
241	WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :',
242	+ LUXLEV,' P=', NSKIP+24
243	ELSE
244	NSKIP = LUXLEV - 24
245	WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV
246	ENDIF
247	IN24 = 0
248	IF (INS .LT. 0) WRITE (6,'(A)')
249	+ ' Illegal initialization by RLUXGO, negative input seed'
250	IF (INS .GT. 0) THEN
251	JSEED = INS
252	WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',
253	+ JSEED, K1,K2
254	ELSE
255	JSEED = JSDFLT_set
256	WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
257	ENDIF
258	INSEED = JSEED
259	NOTYET = .FALSE.
260	TWOM24 = 1.
261	DO 325 I= 1, 24
262	TWOM24 = TWOM24 * 0.5
263	K = JSEED/53668
264	JSEED = 40014(JSEED-K53668) -K*12211
265	IF (JSEED .LT. 0) JSEED = JSEED+ICONS
266	ISEEDS(I) = MOD(JSEED,ITWO24)
267	325 CONTINUE
268	TWOM12 = TWOM24 * 4096.
269	DO 350 I= 1,24
270	SEEDS(I) = REAL(ISEEDS(I))*TWOM24
271	NEXT(I) = I-1
272	350 CONTINUE
273	NEXT(1) = 24
274	I24 = 24
275	J24 = 10
276	CARRY = 0.
277	IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
278	C If restarting at a break point, skip K1 + IGIGA*K2
279	C Note that this is the number of numbers delivered to
280	C the user PLUS the number skipped (if luxury .GT. 0).
281	KOUNT = K1
282	MKOUNT = K2
283	IF (K1+K2 .NE. 0) THEN
284	DO 500 IOUTER= 1, K2+1
285	INNER = IGIGA
286	IF (IOUTER .EQ. K2+1) INNER = K1
287	DO 450 ISK= 1, INNER
288	UNI = SEEDS(J24) - SEEDS(I24) - CARRY
289	IF (UNI .LT. 0.) THEN
290	UNI = UNI + 1.0
291	CARRY = TWOM24
292	ELSE
293	CARRY = 0.
294	ENDIF
295	SEEDS(I24) = UNI
296	I24 = NEXT(I24)
297	J24 = NEXT(J24)
298	450 CONTINUE
299	500 CONTINUE
300	C Get the right value of IN24 by direct calculation
301	IN24 = MOD(KOUNT, NSKIP+24)
302	IF (MKOUNT .GT. 0) THEN
303	IZIP = MOD(IGIGA, NSKIP+24)
304	IZIP2 = MKOUNT*IZIP + IN24
305	IN24 = MOD(IZIP2, NSKIP+24)
306	ENDIF
307	C Now IN24 had better be between zero and 23 inclusive
308	IF (IN24 .GT. 23) THEN
309	WRITE (6,'(A/A,3I11,A,I5)')
310	+ ' Error in RESTARTING with RLUXGO:',' The values', INS,
311	+ K1, K2, ' cannot occur at luxury level', LUXLEV
312	IN24 = 0
313	ENDIF
314	ENDIF
315	RETURN
316	END