C Various CERNLIB mathematical routines needed for the optimum interval,
C as modified from 2006_src.tar.gz to not call other cernlib routines.
C See http://cernlib.web.cern.ch/cernlib/download/2006_source/tar/
C
FUNCTION DFREQ(X)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C
C
DIMENSION P1(0:3),Q1(0:3),P2(0:7),Q2(0:7),P3(0:4),Q3(0:4)
PARAMETER(Z1 = 1, HF = Z1/2)
PARAMETER(C1 = 0.56418 95835 47756 29D0)
PARAMETER(W2 = 1.41421 35623 73095 05D0, RW2 = 1/W2)
DATA (P1(I),Q1(I),I=0,3)
+/+2.42667 95523 05317 5D+2, +2.15058 87586 98612 0D+2,
1 +2.19792 61618 29415 2D+1, +9.11649 05404 51490 1D+1,
2 +6.99638 34886 19135 5D+0, +1.50827 97630 40778 7D+1,
3 -3.56098 43701 81538 5D-2, +1/
DATA (P2(I),Q2(I),I=0,7)
+/+3.00459 26102 01616 01D+2, +3.00459 26095 69832 93D+2,
1 +4.51918 95371 18729 42D+2, +7.90950 92532 78980 27D+2,
2 +3.39320 81673 43436 87D+2, +9.31354 09485 06096 21D+2,
3 +1.52989 28504 69404 04D+2, +6.38980 26446 56311 67D+2,
4 +4.31622 27222 05673 53D+1, +2.77585 44474 39876 43D+2,
5 +7.21175 82508 83093 66D+0, +7.70001 52935 22947 30D+1,
6 +5.64195 51747 89739 71D-1, +1.27827 27319 62942 35D+1,
7 -1.36864 85738 27167 07D-7, +1/
DATA (P3(I),Q3(I),I=0,4)
+/-2.99610 70770 35421 74D-3, +1.06209 23052 84679 18D-2,
1 -4.94730 91062 32507 34D-2, +1.91308 92610 78298 41D-1,
2 -2.26956 59353 96869 30D-1, +1.05167 51070 67932 07D+0,
3 -2.78661 30860 96477 88D-1, +1.98733 20181 71352 56D+0,
4 -2.23192 45973 41846 86D-2, +1/
V=RW2*DABS(X)
IF(V .LT. HF) THEN
Y=V**2
AP=P1(3)
AQ=Q1(3)
DO 1 I = 2,0,-1
AP=P1(I)+Y*AP
1 AQ=Q1(I)+Y*AQ
H=V*AP/AQ
HC=1-H
ELSEIF(V .LT. 4) THEN
AP=P2(7)
AQ=Q2(7)
DO 2 I = 6,0,-1
AP=P2(I)+V*AP
2 AQ=Q2(I)+V*AQ
HC=DEXP(-V**2)*AP/AQ
H=1-HC
ELSE
Y=1/V**2
AP=P3(4)
AQ=Q3(4)
DO 3 I = 3,0,-1
AP=P3(I)+Y*AP
3 AQ=Q3(I)+Y*AQ
HC=DEXP(-V**2)*(C1+Y*AP/AQ)/V
H=1-HC
ENDIF
IF(X .GT. 0) THEN
DFREQ=HF+HF*H
ELSE
DFREQ=HF*HC
ENDIF
RETURN
END
FUNCTION DGAUSN(P)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C Computes a "Normal Deviate"
C Based on G.W. Hill & A.W. Davis, Algorithm 442 Normal Deviate
C Collected Algorithms from CACM
C
PARAMETER (C = 2.50662 82746 31000 50D0)
PARAMETER (Z1 = 1, HF = Z1/2, C1 = 3*Z1/4, C2 = 7*Z1/8, C3 = Z1/3)
IF(P .LE. 0 .OR. P .GE. 1) THEN
H=0
WRITE(6,101) P
ELSEIF(P .EQ. HF) THEN
H=0
ELSE
X=P
IF(P .GT. HF) X=1-P
X=DSQRT(-2.D0*DLOG(X))
X=X-((7.47395D0*X+494.877D0)*X+1637.720D0)/
1 (((X+117.9407D0)*X+908.401D0)*X+659.935D0)
IF(P .LT. HF) X=-X
S=X**2
Z=C*(P-DFREQ(X))*DEXP(HF*S)
H=(((((C1*S+C2)*Z+X)*X+HF)*C3*Z+HF*X)*Z+1)*Z+X
ENDIF
DGAUSN=H
RETURN
101 FORMAT('ARGUMENT P =',1P,D15.8,' NOT IN RANGE')
END
Real*8 FUNCTION GAMDIS(X,A)
IMPLICIT REAL*8 (A-H,O-Z)
C Calculates the gamma distribution function
C
C G(x,a) = (1/gamma(a)) * int(0,x)[exp(-t) * t**(a-1)] dt.
C
C Based on
C W. Gautschi, ALGORITHM 542 Incomplete Gamma Functions,
C ACM Trans. Math. Software 5 (1979) 482-489.
DIMENSION C(14)
PARAMETER (EPS = 1E-5, EPS1 = 5E-7)
PARAMETER (ALH = -0.69314 72)
PARAMETER (Z1 = 1, HALF = Z1/2, QUAR = Z1/4)
PARAMETER (C1 = 3*Z1/2, KMAX = 300)
DATA C
1/ 0.5772157,-0.6558781,-0.0420026, 0.1665386,-0.0421977,
2 -0.0096220, 0.0072189,-0.0011652,-0.0002152, 0.0001281,
3 -0.0000201,-0.0000013, 0.0000011,-0.0000002/
HST=0
IF(X .EQ. 0) GO TO 99
IF(X .LT. 0 .OR. A .LE. 0) THEN
WRITE(6,101) X,A
GO TO 99
ELSE
ALX=DLOG(X)
ENDIF
IF(X .LT. QUAR) THEN
ALFA=ALH/ALX
ELSE
ALFA=X+QUAR
ENDIF
IF(A .GT. ALFA) THEN
TERM=1
SUM=1
DO 1 K = 1,KMAX
TERM=X*TERM/(A+K)
SUM=SUM+TERM
IF(DABS(TERM) .LE. EPS*SUM) GO TO 2
1 CONTINUE
GO TO 98
2 HST=SUM*DEXP(A*ALX-X-ALOGAM(1+A))
ELSEIF(X .GT. C1) THEN
P=0
S=1-A
Q=(X+S)*(X-1-A)
R=4*(X+S)
TERM=1
SUM=1
RHO=0
DO 3 K = 2,KMAX
P=P+S
Q=Q+R
R=R+8
S=S+2
T=P*(1+RHO)
RHO=T/(Q-T)
TERM=RHO*TERM
SUM=SUM+TERM
IF(DABS(TERM) .LE. EPS*SUM) GO TO 4
3 CONTINUE
GO TO 98
4 HST=1-(A*SUM/(X+1-A))*DEXP(A*ALX-X-ALOGAM(1+A))
ELSE
IF(A .LT. HALF) THEN
SUM=C(14)
DO 12 K = 13,1,-1
12 SUM=A*SUM+C(K)
GA=-SUM/(1+A*SUM)
Y=A*ALX
IF(DABS(Y) .GE. 1) THEN
U=GA-(DEXP(Y)-1)/A
ELSE
SUM=1
TERM=1
DO 7 K = 2,KMAX
TERM=Y*TERM/K
SUM=SUM+TERM
IF(DABS(TERM) .LE. EPS1*SUM) GO TO 8
7 CONTINUE
GO TO 98
8 U=GA-SUM*ALX
ENDIF
ELSE
U=GAMMA(A)-DEXP(A*ALX)/A
ENDIF
P=A*X
Q=A+1
R=A+3
TERM=1
SUM=1
DO 9 K = 2,KMAX
P=P+X
Q=Q+R
R=R+2
TERM=-P*TERM/Q
SUM=SUM+TERM
IF(DABS(TERM) .LE. EPS1*SUM) GO TO 10
9 CONTINUE
GO TO 98
10 HST=1-A*(U+SUM*DEXP((1+A)*ALX)/(1+A))/GAMMA(1+A)
ENDIF
99 GAMDIS=HST
RETURN
98 WRITE(6,102) X,A
GO TO 99
101 FORMAT('ILLEGAL ARGUMENT(S) X = ',E15.8,' A = ',E15.8)
102 FORMAT('PROBLEMS WITH CONVERGENCE, X = ',E15.8,' A = ',E15.8)
END
SUBROUTINE RZERO(A,B,X0,R,EPS,MXF,F)
IMPLICIT REAL*8 (A-H,O-Z)
EXTERNAL F
PARAMETER (ONE = 1, HALF = ONE/2)
XA=MIN(A,B)
XB=MAX(A,B)
FA=F(XA,1)
FB=F(XB,2)
IF(FA*FB .GT. 0) GO TO 5
MC=0
1 X0=HALF*(XA+XB)
R=X0-XA
EE=EPS*(DABS(X0)+1)
IF(R .LE. EE) GO TO 4
F1=FA
X1=XA
F2=FB
X2=XB
2 FX=F(X0,2)
MC=MC+1
IF(MC .GT. MXF) GO TO 6
IF(FX*FA .GT. 0) THEN
XA=X0
FA=FX
ELSE
XB=X0
FB=FX
END IF
3 U1=F1-F2
U2=X1-X2
U3=F2-FX
U4=X2-X0
IF(U2 .EQ. 0 .OR. U4 .EQ. 0) GO TO 1
F3=FX
X3=X0
U1=U1/U2
U2=U3/U4
CA=U1-U2
CB=(X1+X2)*U2-(X2+X0)*U1
CC=(X1-X0)*F1-X1*(CA*X1+CB)
IF(CA .EQ. 0) THEN
IF(CB .EQ. 0) GO TO 1
X0=-CC/CB
ELSE
U3=CB/(2*CA)
U4=U3*U3-CC/CA
IF(U4 .LT. 0) GO TO 1
X0=-U3+SIGN(DSQRT(U4),X0+U3)
END IF
IF(X0 .LT. XA .OR. X0 .GT. XB) GO TO 1
R=MIN(DABS(X0-X3),DABS(X0-X2))
EE=EPS*(DABS(X0)+1)
IF(R .GT. EE) THEN
F1=F2
X1=X2
F2=F3
X2=X3
GO TO 2
END IF
FX=F(X0,2)
IF(FX .EQ. 0) GO TO 4
IF(FX*FA .LT. 0) THEN
XX=X0-EE
IF(XX .LE. XA) GO TO 4
FF=F(XX,2)
FB=FF
XB=XX
ELSE
XX=X0+EE
IF(XX .GE. XB) GO TO 4
FF=F(XX,2)
FA=FF
XA=XX
END IF
IF(FX*FF .GT. 0) THEN
MC=MC+2
IF(MC .GT. MXF) GO TO 6
F1=F3
X1=X3
F2=FX
X2=X0
X0=XX
FX=FF
GO TO 3
END IF
4 R=EE
FF=F(X0,3)
RETURN
5 Continue
WRITE(*,100)
R=-2*(XB-XA)
X0=0
RETURN
6 Continue
WRITE(*,101)
R=-HALF*DABS(XB-XA)
X0=0
RETURN
100 FORMAT(1X,'***** Modified RZERO ... F(A) AND F(B)',
1 ' HAVE THE SAME SIGN')
101 FORMAT(1X,'***** Modified RZERO ... TOO MANY FUNCTION CALLS')
END
Real*8 FUNCTION ALGAMA(X)
C
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
Dimension P1(7),Q1(7),P2(7),Q2(7),P3(7),Q3(7),C(5)
PARAMETER (Z1 = 1, HF = Z1/2, HF1 = 1+HF)
DATA P1
1/+3.84287 36567 45991D+0, +5.27068 93753 00983D+1,
2 +5.55840 45723 51531D+1, -2.15135 13573 72570D+2,
3 -2.45872 61722 29242D+2, -5.75008 93603 04123D+1,
4 -2.33590 98949 51284D+0/
DATA Q1
1/+1.00000 00000 00000D+0, +3.37330 47907 07074D+1,
2 +1.93877 84034 37713D+2, +3.08829 54973 42428D+2,
3 +1.50068 39064 89095D+2, +2.01068 51344 33395D+1,
4 +4.57174 20282 50299D-1/
DATA P2
1/+4.87402 01396 83863 6D+0, +2.48845 25168 57407 6D+2,
2 +2.17973 66058 89591 5D+3, +3.79751 24011 52511 8D+3,
3 -1.97780 70769 84164 6D+3, -3.69298 34005 59128 2D+3,
4 -5.60177 73537 80387 7D+2/
DATA Q2
1/+1.00000 00000 00000 0D+0, +9.50999 17418 20893 8D+1,
2 +1.56120 45277 92863 5D+3, +7.23400 87928 94807 1D+3,
3 +1.04595 76594 05895 9D+4, +4.16994 15153 20023 1D+3,
4 +2.76785 83623 80410 1D+2/
DATA P3
1/-6.88062 40094 59425D+3, -4.30699 69819 57098D+5,
2 -4.75045 94653 43956D+6, -2.94234 45930 32234D+6,
3 +3.63218 04931 54257D+7, -3.35677 82814 54576D+6,
4 -2.48043 69488 28593D+7/
DATA Q3
1/+1.00000 00000 00000D+0, -1.42168 29839 65146D+3,
2 -1.55528 90280 85353D+5, -3.41525 17108 01107D+6,
3 -2.09696 23255 80444D+7, -3.45441 75093 34395D+7,
4 -9.16055 82863 71317D+6/
DATA C
1/ 1.12249 21356 561D-1, 7.95916 92961 204D-2,
1 -1.70877 94611 020D-3, 9.18938 53320 467D-1,
2 1.34699 05627 879D+0/
ENTRY ALOGAM(X)
IF(X .LE. 0) THEN
H=0
WRITE(6,101) X
ELSE IF(X .EQ. 1 .OR. X .EQ. 2) THEN
H=0
ELSE IF(X .LE. HF) THEN
Y=X+1
AP=P1(1)
AQ=Q1(1)
DO 2 I = 2,7
AP=P1(I)+Y*AP
2 AQ=Q1(I)+Y*AQ
H=-DLOG(X)+X*AP/AQ
ELSE IF(X .LE. HF1) THEN
AP=P1(1)
AQ=Q1(1)
DO 3 I = 2,7
AP=P1(I)+X*AP
3 AQ=Q1(I)+X*AQ
H=(X-1)*AP/AQ
ELSE IF(X .LE. 4) THEN
AP=P2(1)
AQ=Q2(1)
DO 4 I = 2,7
AP=P2(I)+X*AP
4 AQ=Q2(I)+X*AQ
H=(X-2)*AP/AQ
ELSE IF(X .LE. 12) THEN
AP=P3(1)
AQ=Q3(1)
DO 5 I = 2,7
AP=P3(I)+X*AP
5 AQ=Q3(I)+X*AQ
H=AP/AQ
ELSE
Y=1/X**2
H=(X-HF)*DLOG(X)-X+C(4)+(C(1)+Y*(C(2)+Y*C(3)))/
1 ((C(5)+Y)*X)
ENDIF
ALGAMA=H
RETURN
101 FORMAT('NON-POSITIVE ARGUMENT X = ',1P,E15.6)
END
Real*8 FUNCTION GAMMA(X)
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION C(0:15)
DATA C( 0) /3.65738 77250 83382 44D0/
DATA C( 1) /1.95754 34566 61268 27D0/
DATA C( 2) /0.33829 71138 26160 39D0/
DATA C( 3) /0.04208 95127 65575 49D0/
DATA C( 4) /0.00428 76504 82129 09D0/
DATA C( 5) /0.00036 52121 69294 62D0/
DATA C( 6) /0.00002 74006 42226 42D0/
DATA C( 7) /0.00000 18124 02333 65D0/
DATA C( 8) /0.00000 01096 57758 66D0/
DATA C( 9) /0.00000 00059 87184 05D0/
DATA C(10) /0.00000 00003 07690 81D0/
DATA C(11) /0.00000 00000 14317 93D0/
DATA C(12) /0.00000 00000 00651 09D0/
DATA C(13) /0.00000 00000 00025 96D0/
DATA C(14) /0.00000 00000 00001 11D0/
DATA C(15) /0.00000 00000 00000 04D0/
U=X
IF(U .LE. 0) THEN
WRITE(6,101) U
H=0
B0=0
B2=0
F=0
GO TO 9
ENDIF
8 F=1
IF(U .LT. 3) THEN
DO 1 I = 1,INT(4-U)
F=F/U
1 U=U+1
ELSE
DO 2 I = 1,INT(U-3)
U=U-1
2 F=F*U
END IF
H=U+U-7
ALFA=H+H
B1=0
B2=0
DO 3 I = 15,0,-1
B0=C(I)+ALFA*B1-B2
B2=B1
3 B1=B0
9 GAMMA =F*(B0-H*B2)
RETURN
101 FORMAT('ARGUMENT IS NEGATIVE = ',1P,E15.1)
END
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