c-------------------------------------------------------------------- 
c                    Test runlux routines
c-------------------------------------------------------------------- 
      PARAMETER (Ntot =20 ) ! total number of rand numbs
      DIMENSION RanNum(5)   ,RR(Ntot)
      real*8 sig, xmean
      iseed = 9121071
      CALL rluxgo(2,iseed,0,0)
      icurr = 0
 10   call ranlux(RanNum,5)
c      icurr=icurr +1
         do i=1,5
            icurr=icurr +1
            RR(icurr) = RanNum(i)
         EndDo 
         If(icurr.le.Ntot-5)goto 10

         write (*,*)  ' done'
c         open(1,file='set2.dat')
c         write (1,*) RR
         sig =0.
         xmean =0.
         i1 = 0
         im1 = 0
         i15 = 0
         im15 = 0
         i2 = 0
         im2 = 0
         i25 = 0
         im25 = 0
         i3 = 0
         im3 = 0
         i4 = 0
         im4 = 0
         Ntry =1e8
         s1  = 0.15866
         s15 = 6.681e-2
         s2   = 2.275e-2
         s25 = 6.21e-3
         s3  = 1.35e-3
         s4  = 3.17e-5
         Do i=1,Ntry
c           xx =  GAUSS(iseed)
            xx =  GAUSS3()
            sig = sig +xx**2
            xmean = xmean + xx
            if(xx.gt.1.d+0)i1 =i1+1
            if(xx.lt.-1.d+0)im1 =im1+1
            if(xx.gt.1.5d+0)i15 =i15+1
            if(xx.lt.-1.5d+0)im15 =im15+1
            if(xx.gt.2.d+0)i2 =i2+1
            if(xx.lt.-2.d+0)im2 =im2+1
            if(xx.lt.-2.5d+0)im25 =im25+1
            if(xx.gt.2.5d+0)i25 =i25+1
            if(xx.gt.3.d+0)i3 =i3+1
            if(xx.lt.-3.d+0)im3 =im3+1
            if(xx.gt.4.d+0)i4 =i4+1
            if(xx.lt.-4.d+0)im4 =im4+1
         EndDo
         write (*,*)  ' N=    ',Ntry,' GAUSS3+ranlux ------'
         at =float (Ntry)
         write (*,'("  sigma=",g12.5," mean=",g12.5)') sig/at,xmean/at
         write (*,*)  ' sigma Frac(>sig) Frac(<-sig)   True',
     &                      '       n(>)    n(<)'
         write (*,'(f6.2,3g11.4,2i9)') 1.0,i1/at, im1/at,s1,i1,im1
         write (*,'(f6.2,3g11.4,2i9)') 1.5,i15/at, im15/at,s15,i15,im15
         write (*,'(f6.2,3g11.4,2i9)') 2.0,i2/at, im2/at,s2,i2,im2
         write (*,'(f6.2,3g11.4,2i9)') 2.5,i25/at, im25/at,s25,i25,im25
         write (*,'(f6.2,3g11.4,2i9)') 3.0,i3/at, im3/at,s3,i3,im3
         write (*,'(f6.2,3g11.4,2i9)') 4.0,i4/at, im4/at,s4,i4,im4

         Stop
         End
c-------------------------------------------------------------------- 
            
!   LUXURY LEVELS.
!  level 0  (p=24): equivalent to the original RCARRY of Marsaglia
!           and Zaman, very long period, but fails many tests.
!  level 1  (p=48): considerable improvement in quality over level 0,
!           now passes the gap test, but still fails spectral test.
!  level 2  (p=97): passes all known tests, but theoretically still
!           defective.
!  level 3  (p=223): DEFAULT VALUE.  Any theoretically possible
!           correlations have very small chance of being observed.
!  level 4  (p=389): highest possible luxury, all 24 bits chaotic.
c---------------------------------------------
      SUBROUTINE ranlux(rvec, lenv)
c                    returns a vector RVEC of LEN   
c                   32-bit random floating point numbers between 
c                   zero (not included) and one (also not incl.).
c           Next call to ranlux gives next LEN of random numbers
c---------------------------------------------
      INTEGER  lenv
      REAL        rvec(lenv)
      INTEGER  iseeds(24)
      INCLUDE 'luxury.h'
      DATA ndskip/ 0, 24, 73, 199, 365 /
      DATA i24,j24,luxlev/24,10,lxdflt/
      DATA notyet/.true./
      DATA in24,kount,mkount,carry/0,0,0,0./


!  NOTYET is .TRUE. if no initialization has been performed yet.
!              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 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)
        END DO
        twom12 = twom24 * 4096.
        DO i = 1, 24
          seeds(i) = REAL(iseeds(i)) * twom24
          next(i) = i - 1
        END DO
        next(1) = 24
        i24 = 24
        j24 = 10
        carry = 0.
        IF (seeds(24).EQ.0.) carry = twom24
      END IF

!          The Generator proper: "Subtract-with-borrow",
!          as proposed by Marsaglia and Zaman,
!          Florida State University, March, 1989

      DO ivec = 1, lenv
        uni = seeds(j24) - seeds(i24) - carry
        IF (uni.LT.0.) THEN
          uni = uni + 1.0
          carry = twom24
        ELSE
          carry = 0.
        END IF
        seeds(i24) = uni
        i24 = next(i24)
        j24 = next(j24)
        rvec(ivec) = uni
      !  small numbers (with less than 12 "significant" bits) are "padded".
        IF (uni.LT.twom12) THEN
          rvec(ivec) = rvec(ivec) + twom24 * seeds(j24)
      !        and zero is forbidden in case someone takes a logarithm
          IF (rvec(ivec).EQ.0.) rvec(ivec) = twom24 * twom24
        END IF
      !        Skipping to luxury.  As proposed by Martin Luscher.
        in24 = in24 + 1
        IF (in24.EQ.24) THEN
          in24 = 0
          kount = kount + nskip
          DO isk = 1, nskip
            uni = seeds(j24) - seeds(i24) - carry
            IF (uni.LT.0.) THEN
              uni = uni + 1.0
              carry = twom24
            ELSE
              carry = 0.
            END IF
            seeds(i24) = uni
            i24 = next(i24)
            j24 = next(j24)
          END DO
        END IF
      END DO
      kount = kount + lenv
      IF (kount.GE.igiga) THEN
        mkount = mkount + 1
        kount = kount - igiga
      END IF
      RETURN

      END ! SUBROUTINE ranlux

c---------------------------------------------
c                    Subroutine to initialize from one or three integers
      SUBROUTINE rluxgo(lux, ins, k1, k2)
c                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
c---------------------------------------------
      INCLUDE 'luxury.h'
      INTEGER  iseeds(24)
      
      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 ilx = 0, maxlev
          IF (lux.EQ.ndskip(ilx)+24) luxlev = ilx
        END DO
      END IF
      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
      END IF
      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'
      END IF
      inseed = jseed
      notyet = .false.
      twom24 = 1.
      DO 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)
      END DO
      twom12 = twom24 * 4096.
      DO i = 1, 24
        seeds(i) = REAL(iseeds(i)) * twom24
        next(i) = i - 1
      END DO
      next(1) = 24
      i24 = 24
      j24 = 10
      carry = 0.
      IF (seeds(24).EQ.0.) carry = twom24
      !        If restarting at a break point, skip K1 + IGIGA*K2
      !        Note that this is the number of numbers delivered to
      !        the user PLUS the number skipped (if luxury .GT. 0).
      kount = k1
      mkount = k2
      IF (k1+k2.NE.0) THEN
        DO iouter = 1, k2 + 1
          inner = igiga
          IF (iouter.EQ.k2+1) inner = k1
          DO isk = 1, inner
            uni = seeds(j24) - seeds(i24) - carry
            IF (uni.LT.0.) THEN
              uni = uni + 1.0
              carry = twom24
            ELSE
              carry = 0.
            END IF
            seeds(i24) = uni
            i24 = next(i24)
            j24 = next(j24)
          END DO
        END DO
      !         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)
        END IF
      !       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
        END IF
      END IF
      RETURN
      
      END ! SUBROUTINE rluxgo

C--------------------------------------
C		normal random numbers
      FUNCTION GAUSS1()
C                          Uses ranlux for homogenous rand numbers
C                          sum of 5 random numbers + corrections
C                          It is quite slow: needs 5 calls to ranlux  
C         N =1e8
C  sigma Frac(>sig) Frac(<-sig)   True       n(>)    n(<)
c  1.00 0.1580     0.1580     0.1587          15796874 15798080
c 1.50 0.6674E-01 0.6676E-01 0.6681E-01  6673965  6675533
c  2.00 0.2279E-01 0.2281E-01 0.2275E-01  2278910  2281013
c  2.50 0.6106E-02 0.6108E-02 0.6210E-02   610575   610813
c  3.00 0.1153E-02 0.1151E-02 0.1350E-02   115341   115052
c  4.00 0.2230E-05 0.2070E-05 0.3170E-04      223      207
C--------------------------------------
      DIMENSION RanNum(5) 
      X=0.
      CALL ranlux(RanNum,5)
      DO  I=1,5
         X=X+RanNum(I)
      EndDo
      X2   =1.5491933*(X-2.5)
      GAUSS1=X2  *(1.-0.012*(3.-X2*X2))
      RETURN
      END
C--------------------------------------
C		normal random numbers
      FUNCTION GAUSS2(M)
C                        Uses ranndd for homogenous rand numbers
C                          sum of 12 random numbers + corrections
C                          It is  slow: very (10 times than GAUSS3)
C                          Quality is very good: upto 4 sigma
c  sigma Frac(>sig) Frac(<-sig)   True       n(>)    n(<)
c  1.00 0.1586     0.1586     0.1587     15857116 15857711
c  1.50 0.6678E-01 0.6680E-01 0.6681E-01  6678110  6679725
c  2.00 0.2277E-01 0.2275E-01 0.2275E-01  2277430  2274868
c  2.50 0.6197E-02 0.6189E-02 0.6210E-02   619732   618895
c  3.00 0.1327E-02 0.1319E-02 0.1350E-02   132713   131929
c  4.00 0.2664E-04 0.2727E-04 0.3170E-04     2664     2727
C                   
c-------------------------------------------------------------------- 
      X=0.
      DO  I=1,12
         X=X+Randd(M)
      EndDo
      X2 = X-6.0
      GAUSS2   =X2 *(1.-0.0045*(3.-X2*X2))
      RETURN
      END
C--------------------------------------
C		normal random numbers
      FUNCTION GAUSS3()
C                        Uses ranlux for homogenous rand numbers   
C                        Uses Box-Muller inversion for gaussian
C                 Excellent quality and speed
c  N=       100000000  GAUSS3+ranlux ------
c  sigma= 0.99997     mean= 0.18692E-03
c  sigma Frac(>sig) Frac(<-sig)   True       n(>)    n(<)
c  1.00 0.1587     0.1586     0.1587     15869095 15857724
c  1.50 0.6682E-01 0.6680E-01 0.6681E-01  6681652  6679768
c  2.00 0.2277E-01 0.2273E-01 0.2275E-01  2276651  2273197
c  2.50 0.6222E-02 0.6206E-02 0.6210E-02   622190   620610
c  3.00 0.1356E-02 0.1347E-02 0.1350E-02   135628   134674
c  4.00 0.3242E-04 0.3145E-04 0.3170E-04     3242     3145
c
C--------------------------------------
      DIMENSION RanNum(2) 
      DATA iFlag/0/
      SAVE gSet
      If(iFlag.eq.0)Then
 1       CALL ranlux(RanNum,2)
         x1 = 2.*RanNum(1)-1.
         x2 = 2.*RanNum(2)-1.
         R =x1**2+x2**2
         If(R.ge.1.)GoTo 1
         F =sqrt(-2.*LOG(R)/R)
         gSet    = x1*F
         GAUSS3 = x2*F
         iFlag = 1
      Else
         GAUSS3 = gSet
         iFlag = 0
      EndIf 
      RETURN
      END

C-------------------------------------------               random number generator
      FUNCTION RANDd(M)
C------------------------------------------------
      DATA LC,AM,KI,K1,K2,K3,K4,L1,L2,L3,L4
     +	/453815927,2147483648.,2147483647,536870912,131072,256,
     +	 16777216,4,16384,8388608,128/
      ML=M/K1*K1
      M1=(M-ML)*L1
      ML=M/K2*K2
      M2=(M-ML)*L2
      ML=M/K3*K3
      M3=(M-ML)*L3
      ML=M/K4*K4
      M4=(M-ML)*L4
      M5=KI-M
      IF(M1.GE.M5)M1=M1-KI-1
      ML=M+M1
      M5=KI-ML
      IF(M2.GE.M5)M2=M2-KI-1
      ML=ML+M2
      M5=KI-ML
      IF(M3.GE.M5)M3=M3-KI-1
      ML=ML+M3
      M5=KI-ML
      IF(M4.GE.M5)M4=M4-KI-1
      ML=ML+M4
      M5=KI-ML
      IF(LC.GE.M5)ML=ML-KI-1
      M=ML+LC
      RANDd=M/AM
      RETURN
      END
C--------------------------------------
C		normal random numbers
      FUNCTION GAUSS(M)
C                    Uses RANDd() + 5 random numbers with correction
C                    Old code used in many simulations
C   N=1e8  sigma= 0.99556     mean= 0.10542E-03
c  sigma Frac(>sig) Frac(<-sig)   True       n(>)    n(<)
c  1.00 0.1589     0.1589     0.1587     15894145 15889594
c  1.50 0.6704E-01 0.6702E-01 0.6681E-01  6703736  6702430
c  2.00 0.2258E-01 0.2256E-01 0.2275E-01  2258038  2256030
c  2.50 0.5841E-02 0.5843E-02 0.6210E-02   584148   584264
c  3.00 0.1031E-02 0.1028E-02 0.1350E-02   103141   102783
c  4.00 0.5900E-06 0.6100E-06 0.3170E-04       59       61
C
C--------------------------------------
      X=0.
      DO  I=1,5
         X=X+RANDd(M)
      EndDo
      X2   =1.5491933*(X-2.5)
      GAUSS=X2*(1.-0.01*(3.-X2*X2))
      RETURN
      END