!     Last change:  DOS   3 Aug 2000    5:59 pm
!                                       *** copyright 2000 ***
! *** filename fstmed.f95               *** John F. Monahan **
!                                       **********************
      program pfstmed
!  TEST OF FAST MEDIAN ALGORITHM FSTMED
      implicit none
      real, dimension(30000) :: x
      real xmed,xmeds
      integer i,j,k,n,ka,kb,ke,iflag
      real ran,rang,fstmed    ! rang just a second copy of ranbf
!
  20  FORMAT(//' Test Problems with Sample Size',i8)
  21  FORMAT(/' Simple Case -- All Positive, N=',i8,'  OK if zero',i2)
  23  FORMAT(/' Simple Case -- All Negative, N=',i8,'  OK if zero',i2)
  25  FORMAT(/' Simple Case -- Mixed Signs, N=',i8,'  OK if zero',i2)
  27  FORMAT(/' Mixed Signs and Some Ties, N=',i8,'  OK if zero',i2)
!  22  FORMAT(2X,8F9.4)
  28  FORMAT(2X,8F9.4)
  29  FORMAT(' Median via fstmed',f12.6,4x,'via sorting',f12.6)
!  28  FORMAT(2X,12F6.2)
!                   open output file
      open( unit=6, file='fstmed.out' )
!                   initialize uniform generators
      xmed = ran(5151917)
      xmed = rang(5151917)
!                   let's check out these hundred and twenty cases
      do i = 1,120
      n = i 
      if( i .gt. 10 ) n = i*i
      write(6,20) n
!                   first all positive
      do j = 1,n
      x(j) = 5.*rang(j)
      end do  ! loop on j
      if( n .lt. 20 ) write(6,28) (x(j),j=1,n)
!                   call median algorithm and write out results
      xmed = fstmed(x,n)
!                   test
      ka = 0
      kb = 0
      ke = 0
      do j = 1,n
      if( x(j) .gt. xmed ) ka = ka + 1
      if( x(j) .eq. xmed ) ke = ke + 1
      if( x(j) .lt. xmed ) kb = kb + 1
      end do  ! loop on j
      iflag = 1
      if( (ka+ke .ge. kb) .and. (kb+ke .ge. ka) ) iflag = 0
!                   check against sorting route
      call qsort(x,n)
      xmeds = (x( (n+1)/2 ) + x( (n+2)/2 ) )/2.
      write(6,29) xmed,xmeds
      if( abs(xmed - xmeds) .gt. 1.e-4 ) iflag = 1
      write(6,21) n,iflag
      if( n .lt. 20 ) write(6,28) (x(j),j=1,n)
!                   now all negative
      do j = 1,n
      x(j) = -100.*rang(j)
      end do  ! loop on j
      if( n .lt. 20 ) write(6,28) (x(j),j=1,n)
!                   call median algorithm and write out results
      xmed = fstmed(x,n)
!                   test
      ka = 0
      kb = 0
      ke = 0
      do j = 1,n
      if( x(j) .gt. xmed ) ka = ka + 1
      if( x(j) .eq. xmed ) ke = ke + 1
      if( x(j) .lt. xmed ) kb = kb + 1
      end do  ! loop on j
      iflag = 1
      if( (ka+ke .ge. kb) .and. (kb+ke .ge. ka) ) iflag = 0
!                   check against sorting route
      call qsort(x,n)
      xmeds = (x( (n+1)/2 ) + x( (n+2)/2 ) )/2.
      write(6,29) xmed,xmeds
      if( abs(xmed - xmeds) .gt. 1.e-4 ) iflag = 1
      write(6,23) n,iflag
      if( n .lt. 20 ) write(6,28) (x(j),j=1,n)
!                   now mixed positive and negative and smaller
      do j = 1,n
      x(j) = (rang(j)-.3)/100.
      end do  ! loop on j
      if( n .lt. 20 ) write(6,28) (x(j),j=1,n)
!                   call median algorithm and write out results
      xmed = fstmed(x,n)
!                   test
      ka = 0
      kb = 0
      ke = 0
      do j = 1,n
      if( x(j) .gt. xmed ) ka = ka + 1
      if( x(j) .eq. xmed ) ke = ke + 1
      if( x(j) .lt. xmed ) kb = kb + 1
      end do  ! loop on j
      iflag = 1
      if( (ka+ke .ge. kb) .and. (kb+ke .ge. ka) ) iflag = 0
!                   check against sorting route
      call qsort(x,n)
      xmeds = (x( (n+1)/2 ) + x( (n+2)/2 ) )/2.
      write(6,29) xmed,xmeds
      if( abs(xmed - xmeds) .gt. 1.e-4 ) iflag = 1
      write(6,25) n,iflag
      if( n .lt. 20 ) write(6,28) (x(j),j=1,n)
!                   now introduce some ties
      do j = 1,n
      k = int( 8.*(rang(j)-.3) )
      x(j) = real(k)/10.
      end do  ! loop on j
      if( n .lt. 20 ) write(6,28) (x(j),j=1,n)
!                   call median algorithm and write out results
      xmed = fstmed(x,n)
!                   test
      ka = 0
      kb = 0
      ke = 0
      do j = 1,n
      if( x(j) .gt. xmed ) ka = ka + 1
      if( x(j) .eq. xmed ) ke = ke + 1
      if( x(j) .lt. xmed ) kb = kb + 1
      end do  ! loop on j
      iflag = 1
      if( (ka+ke .ge. kb) .and. (kb+ke .ge. ka) ) iflag = 0
!                   check against sorting route
      call qsort(x,n)
      xmeds = (x( (n+1)/2 ) + x( (n+2)/2 ) )/2.
      write(6,29) xmed,xmeds
      if( abs(xmed - xmeds) .gt. 1.e-4 ) iflag = 1
      write(6,27) n,iflag
      if( n .lt. 20 ) write(6,28) (x(j),j=1,n)
!                   done with this sample size
      end do  ! loop on i (n)
      stop
      end program pfstmed
      REAL FUNCTION FSTMED(X,N)
!  FAST ALGORITHM FOR COMPUTING MEDIAN OF SAMPLE OF SIZE N
!
!  ARGUMENTS
!  X   REAL VECTOR OF OBSERVATIONS
!  N   INTEGER NUMBER OF OBSERVATIONS
!
!  *** REQUIRED SUBPROGRAMS ***
!  SUBROUTINE PARTIT(X,LEFT,RGHT,NS,NT)    PARTITIONS SET
!  REAL FUNCTION RAN(IDUM)                 UNIFORM(0,1) RANDOM NUMBERS
!
!  J F MONAHAN 
!     RECODED MARCH, APRIL 2000 FOR FORTRAN 95
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: N
      REAL, DIMENSION(N), INTENT(IN OUT) :: X
      INTEGER J,K,KP,START,SPAN,NS,NT
      REAL RAN
!
!                   FIND MEDIAN -- K AND KP MAY OR MAY NOT BE DIFFERENT
      K = (N + 1) / 2
      KP = (N + 2) / 2
!                   INITIAL VALUES
      START = 1
      SPAN = N
!                   RESTART HERE
      DO            ! *** UNRESTRICTED DO ***
!                   RANDOM PARTITION ELEMENT (UNLESS TOO SHORT)
      NS = START
      IF( SPAN .GT. 2 ) NS = START + INT( RAN(SPAN)*REAL(SPAN) )
!                   PARTITION X FROM START TO START+SPAN-1 WITH X(NS)
      CALL PARTIT(X,START,START+SPAN-1,NS,NT)
!                   NS = NUMBER SMALLER, NT = NUMBER TIED
!
!                   PARTITION MAKES TWO SUBLISTS 
!
!                   ARE WE NEARLY DONE? 
      IF( ( K .GE. START+NS-1 ) .AND. ( KP .LE. START+NS+NT ) ) EXIT
!                   LOOK AMONG SMALLER
      IF( K .LT. START + NS - 1 ) SPAN = NS
!                   LOOK AMONG LARGER
      IF( KP .GT. START + NS + NT ) THEN
      SPAN = SPAN - NS - NT
      START = START + NS + NT
                                    END IF ! ( KP .GT. START + NS + NT )
      END DO  ! *** END UNRESTRICTED DO ***
!                   NEARLY DONE -- MAY HAVE TO FIND MIN OR MAX
      IF( KP .EQ. START+NS+NT ) THEN
!                   FIND MIN OF LARGERS
      FSTMED = X(START+NS+NT)
      SPAN = SPAN - NS - NT
      DO J = 1,SPAN
      FSTMED = MIN( FSTMED, X(START+NS+NT+J-1) )
      END DO  ! LOOP ON J
      IF( K .NE. KP ) FSTMED = ( FSTMED + X(START+NS) ) / 2.
      RETURN
                                END IF ! ( KP .EQ. START+NS+NT )
!
      IF( K .EQ. START+NS-1 ) THEN
!                   FIND MAX OF SMALLERS
      FSTMED = X(START)
      DO J = 1,NS
      FSTMED = MAX( FSTMED, X(START+J-1) )
      END DO  ! LOOP ON J
      IF( K .NE. KP ) FSTMED = ( FSTMED + X(START+NS) ) / 2.
      RETURN
                              END IF ! ( K .EQ. START+NS-1 )
!                   DONE -- MEDIAN IN MIDDLE
      FSTMED = X(START+NS)
      RETURN
      END FUNCTION FSTMED
      SUBROUTINE QSORT(K,N)
!  QUICKSORT 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
!
!  *** REQUIRED SUBPROGRAMS ***
!  SUBROUTINE PARTIT(X,LEFT,RGHT,NS,NT)    PARTITIONS SET
!  REAL FUNCTION RAN(IDUM)                 UNIFORM(0,1) RANDOM NUMBERS
!
!  *** HAS STACK LENGTH OF ONLY 100 *** SHOULD BE PLENTY ***
!  J F MONAHAN
!     RECODED MARCH, APRIL 2000 FOR FORTRAN 95
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: N
      REAL, DIMENSION(N) :: K
      INTEGER, DIMENSION(100) :: STARTS,SPANS
      INTEGER LSTACK,START,SPAN,NS,NT
      REAL RAN
!                   START STACK
      LSTACK = 1
      STARTS(LSTACK) = 1
      SPANS(LSTACK) = N
!                   RESTART HERE
      DO WHILE ( LSTACK .GT. 0 )
      START = STARTS(LSTACK)
      SPAN = SPANS(LSTACK)
!                   RANDOM PARTITION ELEMENT (UNLESS TOO SHORT)
      NS = START
      IF( SPAN .GT. 2 ) NS = START + INT( RAN(LSTACK)*REAL(SPAN) )
!                   PARTITION X FROM START TO START+SPAN-1 WITH X(NS)
      CALL PARTIT(K,START,START+SPAN-1,NS,NT)
      LSTACK = LSTACK - 1
!                   NS = NUMBER SMALLER, NT = NUMBER TIED
!
!                   PARTITION MAKES TWO SUBLISTS, STORE IN STACK 
!
!                   THOSE SMALLER THAN PARTITION ELEMENT X(NS) (IN)
      IF( NS .GT. 1 ) THEN
      LSTACK = LSTACK + 1
      STARTS(LSTACK) = START
      SPANS(LSTACK) = NS
                      END IF ! ( NS .GT. 1 )
!                   THOSE LARGER THAN PARTITION ELEMENT X(NS) (IN)      
      SPAN = SPAN - NS - NT
      IF( SPAN .GT. 1 ) THEN
      LSTACK = LSTACK + 1
      STARTS(LSTACK) = START + NS + NT
      SPANS(LSTACK) = SPAN
                        END IF ! ( SPAN .GT. 1 )
!                   STACK LENGTH OF 100 SHOULD BE GRACIOUS PLENTY
      IF( LSTACK .GT. 100 ) THEN
      WRITE(*) 'STACK OVERFLOW IN QSORT'
      STOP
                            END IF ! ( LSTACK .GT. 100 )
      END DO  ! WHILE ( LSTACK .GT. 0 )
      RETURN
      END SUBROUTINE QSORT
      SUBROUTINE PARTIT(X,LEFT,RGHT,NS,NT)
!  ALGORITHM FOR PARTITIONING LIST OF KEYS K FROM X(LEFT) TO X(RGHT)
!   INTO THREE GROUPS: THOSE < X*, THOSE = X*, THOSE > X*
!
!  ARGUMENTS
!  X      IN    REAL     LIST OF KEYS TO BE PARTITIONED
!         OUT            PARTITIONED LIST
!                        X(LEFT) ... X(LEFT+NS-1) < X*
!                        X(LEFT+NS) ... X(LEFT+NS+NT-1) = X*
!                        X(LEFT+NS+NT) ... X(RGHT) > X*
!
!  LEFT   IN    INTEGER  LEFT ENDPOINT OF LIST
!  RGHT   IN    INTEGER  RIGHT ENDPOINT OF LIST
!  NS     IN    INTEGER  POINTS TO PARTITION ELEMENT X* = X(NS) (IN)
!         OUT            NUMBER IN LIST < X*
!  NT     OUT   INTEGER  NUMBER IN LIST = X*
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: LEFT,RGHT
      INTEGER, INTENT(IN OUT) :: NS
      INTEGER, INTENT(OUT) :: NT
      REAL, DIMENSION(RGHT) :: X
      REAL XCUT,XX
      INTEGER I,J
!                   USE X(NS) TO SPLIT X
      XCUT = X(NS)
!                   DO NOTHING IF THERE'S NOTHING TO DO
      NS = 0
      NT = 1
      IF( RGHT - LEFT .LE. 0 ) RETURN
!                   FIRST PARITITION TO TWO GROUPS: LT AND GE
!
!                   INITIALIZE POINTERS I (LEFT) AND J (RIGHT)
      I = LEFT
      J = RGHT
      DO WHILE( I .LE. J )
      DO WHILE( X(I) .LT. XCUT )
      I = I + 1
      IF( I .GT. J ) EXIT     ! LAST ONE LT IS I-1
      END DO  ! WHILE ( X(I) .LT. XCUT )
      IF( I .GT. J ) EXIT     ! FINISH
!                   HAVE X(I) OUT OF PLACE
      DO WHILE ( X(J) .GE. XCUT )
      J = J - 1
      IF( I .GE. J ) EXIT
      END DO  ! WHILE ( X(J) .GE. XCUT )
      IF( I .GE. J ) EXIT
!                   HAVE X(J) OUT OF PLACE
!                   EXCHANGE I AND J
      XX = X(I)
      X(I) = X(J)
      X(J) = XX
      I = I + 1
      J = J - 1
      END DO  ! FIRST MAIN WHILE
!                   HALFWAY DONE -- X(I-1) IS END OF SMALLER
      NS = I - LEFT
!                   NOW SEPARATE THE TIED FROM THE LARGER
!
!
!                   INITIALIZE POINTERS I (LEFT) AND J (RIGHT)
      J = RGHT
      DO WHILE( I .LE. J )
      DO WHILE( X(I) .EQ. XCUT )
      I = I + 1
      IF( I .GT. J ) EXIT     ! LAST ONE EQ IS I-1
      END DO  ! WHILE ( X(I) .EQ. XCUT )
      IF( I .GT. J ) EXIT     ! FINISH
!                   HAVE X(I) OUT OF PLACE
      DO WHILE( X(J) .GT. XCUT )
      J = J - 1
      IF( I .GE. J ) EXIT
      END DO  ! WHILE ( X(J) .GT. XCUT )
      IF( I .GE. J ) EXIT
!                   HAVE X(J) OUT OF PLACE
!                   EXCHANGE I AND J
      XX = X(I)
      X(I) = X(J)
      X(J) = XX
      I = I + 1
      J = J - 1
      END DO  ! SECOND MAIN WHILE
      NT = I - LEFT - NS
      RETURN
      END SUBROUTINE PARTIT
      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
      REAL FUNCTION RANG(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
      RANG = REAL(IX)*RP
      RETURN
      END FUNCTION RANG
! *** end of filename fstmed.f95        *********************