!     Last change:  DOS   1 Aug 2000    8:07 pm
!                                       *** copyright 2000 ***
! *** filename gevonn.f95               *** John F. Monahan **
!                                       **********************
      PROGRAM PGEVONN
!  TEST OF VON NEUMANN ALGORITHM FOR EXPONENTIAL DISTRIBUTION
!                   let's check out 2**14 = 16384
      IMPLICIT NONE
      INTEGER, PARAMETER :: N = 16384
      INTEGER I
      REAL, DIMENSION(N) :: X
      REAL XS,X2S,FN,PHIXI,DKSP,DKSM,DKS,AD
      REAL RAN,GEVONN
!
  21  FORMAT(' Sample Mean & Variance for',I8,' Exponential Deviates' &
     &  /'  Mean is',F16.8,4x,'Variance is',F16.8)
  22  FORMAT(' Goodness of Fit Statistics: '/4x,'Anderson-Darling is',&
     &F12.6/4x,'Compare to .05, .01 critical values of 2.492, 3.857'//&
     &  4x,'Kolmogorov-Smirnov   D+',F9.6,4x,'  D-',F9.6/             &
     &  4x,'SQRT(N)*MAX(D+,D-) is',F9.4/4x,'Compare to .05, .01',     &
     &  ' critical values of 1.36, 1.63')
!
!                   open output file
      OPEN( UNIT=6, FILE='gevonn.out' )
!                   let's check out 2**14 = 16384
      FN = REAL(N)
!                   initialize uniform generator
      XS = RAN(5151917)
!                   initialize moment sums
      X2S = 0.
      XS = 0.
      DO I = 1,N
      X(I) = GEVONN(I)
      XS = XS + X(I)
      IF( I .GT. 1 ) X2S = X2S + (XS - I*X(I))*(XS - I*X(I))/(I*(I-1))
      END DO  ! LOOP ON I
!                   get sample mean and variance
      XS = XS/FN
      X2S = X2S/(FN-1.)
      WRITE(6,21) N,XS,X2S
!                   sort the data -- get order statistics
      CALL HSORT(X,N)
!                   now compute K-S and A-D statistics
      DKSP = 0.
      DKSM = 0.
      AD = 0.
      DO I = 1,N
      PHIXI = 1. - EXP( - X(I) )
      AD = AD - REAL(2*I-1)*( LOG( PHIXI ) - X(N+1-I) )/FN - 1.
      DKSP = MAX( DKSP, REAL(I)/FN - PHIXI )
      DKSM = MAX( DKSM, PHIXI - REAL(I-1)/FN )
      END DO  ! LOOP ON I
!                   usual standardization of K-S
      DKS = SQRT(FN)*MAX( DKSP, DKSM )
      WRITE(6,22) AD,DKSP,DKSM,DKS
      STOP
      END PROGRAM PGEVONN
      REAL FUNCTION GEVONN(IX)
!  VON NEUMANN ALGORITHM FOR GENERATING EXPONENTIAL RANDOM VARIABLES
!
!  ARGUMENT  IX  IS A DUMMY AND NEVER USED
!
!  JOHN VON NEUMANN (1951) 'VARIOUS TECHNIQUES IN CONNECTION WITH
!   RANDOM DIGITS,' IN MONTE CARLO METHOD, NATIONAL BUREAU OF
!   STANDARDS, APPLIED MATH SERIES 12, PP. 36-38
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: IX
      REAL RAN
      REAL U1,U,US
!                   INITIALIZE TO ZERO
      GEVONN = 0.
!                   RESTART HERE BY GENERATING FIRST UNIFORM
      DO          ! NOTE UNCONDITIONAL DO
      U1 = RAN(1)
!                   SAVE IT
      U = U1
      DO          ! NOTE ANOTHER UNCONDITIONAL DO
!                   GENERATE U(N+1) WITH N ODD
      US = RAN(2)
!                   IF NOT SMALLER THAN PREVIOUS (U) (AND N ODD) * DONE
      IF( U .LT. US ) THEN
!                   ACCEPT ==> ADD SHIFT TO FIRST UNIFORM AND GO
      GEVONN = GEVONN + U1
      RETURN
                      END IF ! ( U .LT. US )
!                   IF SMALLER THAN PREVIOUS, THEN CONTINUE SAMPLING
      U = RAN(3)
!                   IF BIGGER THAN PREVIOUS (US) THEN QUIT SAMPLING
      IF( U .GE. US ) EXIT
      END DO  ! INNER DO
!                   IF NOT SMALLER THAN PREVIOUS (AND N EVEN) * REJECT
      GEVONN = GEVONN + 1.
!                   REJECT ==> INCREMENT SHIFT AND RESTART
      END DO  ! OUTER DO
      END FUNCTION GEVONN
      SUBROUTINE HSORT(K,N)
!  HEAPSORT ALGORITHM FOR SORTING ON VECTOR OF KEYS K OF LENGTH N
!
!  ARGUMENTS
!  K   REAL VECTOR OF KEYS TO BE SORTED
!  N   NUMBER OF ITEMS TO BE SORTED
!
!  TO SORT A PARALLEL VECTOR OF RECORDS, USE HKSORT
!
!  J F MONAHAN (DEC, 1999) FORTRAN 95
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: N
      REAL, DIMENSION(N), INTENT(IN OUT) :: K
      REAL KK
      INTEGER I,L,NCUR
!
!                   DO NOTHING IF THERE'S NOTHING TO DO
      IF( N .LE. 1 ) RETURN
!                   INITIALIZE TO BUILDHEAP PART (LOOP ON L)
      L = N/2 + 1
      NCUR = N
      DO I = L,1,-1
      CALL HEAPIFY(I)
      END DO  ! LOOP ON I
      DO I = 2,N
!                   SWITCH CURRENT LARGEST WITH BOTTOM
      KK = K(1)
      K(1) = K(NCUR)
      K(NCUR) = KK
!                   REHEAP WITH ONE SHORTER
      NCUR = NCUR - 1
      CALL HEAPIFY(1)
      END DO  ! LOOP ON NCUR
      RETURN
      CONTAINS

      SUBROUTINE HEAPIFY(II)
      INTEGER, INTENT(IN) :: II
      INTEGER I,J
      I = II
      DO
      J = 2*I                 ! IS IT A LEAF OR ARE THERE SONS?
      IF( J > NCUR ) EXIT     ! A LEAF
      IF( J < NCUR ) THEN     ! ANOTHER SON OF I
                     IF( K(J+1) > K(J) ) J = J+1  ! LARGER SON IS K
                     END IF   ! A LEAF
      IF( K(J) > K(I) ) THEN  ! EXCHANGE
                        KK = K(J)
                        K(J) = K(I)
                        K(I) = KK
                        I = J
                        ELSE  ! EXIT -- HEAP PROPERTY
                        EXIT
                        END IF
      END DO ! WHILE NOT A LEAF
      END SUBROUTINE HEAPIFY

      END SUBROUTINE HSORT
      REAL FUNCTION RAN(IDUM)
!  PORTABLE IMPLEMENTATION OF UNIFORM PSEUDORANDOM NUMBER GENERATOR
!    LEWIS, GOODMAN, & MILLER MULTIPLICATIVE GENERATOR
!     X(N+1) = MOD( 16807*X(N), 2**31-1 )
!
!  P. BRANTLEY, B.L. FOX, L. SCHRAGE (1983) A GUIDE TO SIMULATION
!   SPRINGER-VERLAG, NEW YORK.  PP. 200-202.
!  UPDATED VERSION OF
!  LINUS SCHRAGE (1979)'A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR' 
!   ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 5, PP. 132-138.
!
!  ARGUMENT 
!   IDUM    INTEGER  FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: IDUM
      REAL, PARAMETER :: RP = 4.656612875E-10     ! 1/P
      INTEGER, PARAMETER :: A = 16807             ! MULTIPLIER
      INTEGER, PARAMETER :: B = 127773            ! B = P / A
      INTEGER, PARAMETER :: C = 2836              ! C = P MOD A
      INTEGER, PARAMETER :: P = 2147483647        ! MODULUS 2**31 - 1
      INTEGER, SAVE :: IX = 0
      INTEGER K1
!
!                   IF NOT FIRST CALL, THEN SKIP SETTING SEED
      IF( IX .EQ. 0) IX = IDUM
!                   WRITE NUMBER AS ALPHA*2**16 + BETA
      K1 = IX / B
      IX = A*( IX - K1*B) - K1*C
!                   ABOVE DOES A*IX MOD B -K1*C
      IF( IX .LT. 0 ) IX = IX + P 
!                   RP BELOW IS RECIPROCAL OF P
      RAN = REAL(IX)*RP
      RETURN
      END FUNCTION RAN
! *** end of filename gevonn.f95        *********************