program main
      implicit none

      integer MAXN
      parameter (MAXN = 100000)

      real*8  value(MAXN),checksum
      integer n,seed,m,i

c User inputs.

 10   write (6, *) 'Enter number of values to sort: '
      read (5, *) n
      if (n .gt. MAXN) goto 10

      write (6, *) 'Enter random number seed: '
      read (5, *) seed

 20   write (6, *) 'Enter which final value to print out: '
      read (5, *) m
      if (m .gt. n) goto 20

c Initialize the values.

      call initialize_values(n,value,seed)

c Sort the values.

      call hpsort(n,value)

c Compute and print output.

      checksum = 0.0
      do i = 1,n
        checksum = checksum + i*value(i)
      enddo
      checksum = checksum/n

      write (6,*) 'Sorted value',m,' =',value(m)
      write (6,*) 'Checksum =',checksum

      end


c Initialize the values randomly.

      subroutine initialize_values(n,value,seed)
      implicit none

      integer   n
      real*8    value(*)
      integer   seed

      integer   i
      real*8    rand1

c generate random values
c random number generator is fully periodic, so will not create duplicates

      do i = 1,n
        value(i) = rand1(seed)
      enddo

      return
      end


c Heap sort from Numerical Recipes.

      subroutine hpsort(n,ra)

      integer n
      real*8 ra(n)

      integer i,ir,j,l
      real*8 rra

      if (n.lt.2) return
      l=n/2+1
      ir=n
10    continue
        if(l.gt.1)then
          l=l-1
          rra=ra(l)
        else
          rra=ra(ir)
          ra(ir)=ra(1)
          ir=ir-1
          if(ir.eq.1)then
            ra(1)=rra
            return
          endif
        endif
        i=l
        j=l+l
20      if(j.le.ir)then
          if(j.lt.ir)then
            if(ra(j).lt.ra(j+1))j=j+1
          endif
          if(rra.lt.ra(j))then
            ra(i)=ra(j)
            i=j
            j=j+j
          else
            j=ir+1
          endif
        goto 20
        endif
        ra(i)=rra
      goto 10

      end


C Park/Miller RNG
 
      double precision function rand1(iseed)
      double precision aa,mm,sseed
      parameter (aa=16807.0D0,mm=2147483647.0D0)
      
      sseed = iseed
      sseed = mod(aa*sseed,mm)
      iseed = sseed
      rand1 = sseed/mm

      return
      end