C      ALGORITHM 780, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 24, NO. 1,     March, 1998, P.  102--106.
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	Drivers/
#	Drivers/driver.f
#	Src/
#	Src/src.f
# This archive created: Thu Oct  1 17:07:59 1998
export PATH; PATH=/bin:$PATH
if test ! -d 'Drivers'
then
	mkdir 'Drivers'
fi
cd 'Drivers'
if test -f 'driver.f'
then
	echo shar: will not over-write existing file "'driver.f'"
else
cat << SHAR_EOF > 'driver.f'
      PROGRAM DRIVER
!
!     The return variable, x, is in COMMON to make certain that no
!     overly-enthuiastic optimizer eliminates the function calls.
!
      common x
!
      n = 1000000
      call msec(mtime0)
      do i=1,n
        x = rexp_lg()
      enddo
      call msec(mtime1)
      dt_lg = 0.001*float(mtime1-mtime0)
      dta_lg = 1.E6*dt_lg/float(n)
      print *, 'Old logarithmic method:'
      print *, '  LG Run time =',dt_lg,' seconds'
      print *, '  LG Average time =',dta_lg
!
      call msec(mtime0)
      do i=1,n
        x = rexpu()
      enddo
      call msec(mtime1)
      dt_eau = 0.001*float(mtime1-mtime0)
      dta_eau = 1.E6*dt_eau/float(n)
      print *, 'Unstructured:'
      print *, '  EA Run time =',dt_eau,' seconds'
      print *, '  EA Average time =',dta_eau
!
      call msec(mtime0)
      do i=1,n
        x = rexps()
      enddo
      call msec(mtime1)
      dt_eas = 0.001*float(mtime1-mtime0)
      dta_eas = 1.E6*dt_eas/float(n)
      print *, 'Structured:'
      print *, '  EA Run time =',dt_eas,' seconds'
      print *, '  EA Average time =',dta_eas
!
      stop
      end
      SUBROUTINE MSEC(M)
      integer m, iv(8)
      call date_and_time(values=iv)     ! F90 intrinsic
      m = iv(8) + 1000*(iv(7) + 60*(iv(6) + 60*iv(5)))
      return
      end
      FUNCTION REXP_LG()                ! Old method for comparison
      call random_number(u)
      rexp_lg = -alog(u)
      return
      end
SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Src'
then
	mkdir 'Src'
fi
cd 'Src'
if test -f 'src.f'
then
	echo shar: will not over-write existing file "'src.f'"
else
cat << SHAR_EOF > 'src.f'
      FUNCTION REXPU()
!
!     Random-number generator for the exponential distribution
!     Algorithm EA from J. H. Ahrens and U. Dieter,
!     Communications of the ACM, 31 (1988) 1330--1337.
!     Coded by K. G. Hamilton, December 1996, with corrections.
!
      data alog2/0.6931 4718 0559 9453/
      data     a/5.7133 6315 2645 4228/
      data     b/3.4142 1356 2373 0950/
      data    c/-1.6734 0532 4028 4925/
      data     p/0.9802 5814 3468 5472/
      data    aa/5.6005 7075 6973 8080/
      data    bb/3.3468 1064 8056 9850/
      data    hh/0.0026 1067 2360 2095/
      data    dd/0.0857 8643 7626 9050/
!
    1 call random_number(u)
!       Comment out the following line if your RNG can never return 0.0
      if (u.le.0) go to 1               ! Zero-protector
      g = c
    2 u = u+u
      if (u.ge.1.0) go to 4
      g = g + alog2
      go to 2
    4 u = u-1.0
      if (u.gt.p) go to 6
      rexpu = g + aa/(bb-u)
      return
    6 call random_number(u)
      y = a/(b-u)
      call random_number(up)
      if ((up*hh+dd)*(b-u)**2 .gt. exp(-(y+c))) go to 6
      rexpu = g+y
      return
      end
      FUNCTION REXPS()
!
!     Random-number generator for the exponential distribution
!     Algorithm EA from J. H. Ahrens and U. Dieter,
!     Communications of the ACM, 31 (1988) 1330--1337.
!     Coded by K. G. Hamilton, December 1996, with corrections.
!
      data alog2/0.6931 4718 0559 9453/
      data     a/5.7133 6315 2645 4228/
      data     b/3.4142 1356 2373 0950/
      data    c/-1.6734 0532 4028 4925/
      data     p/0.9802 5814 3468 5472/
      data    aa/5.6005 7075 6973 8080/
      data    bb/3.3468 1064 8056 9850/
      data    hh/0.0026 1067 2360 2095/
      data    dd/0.0857 8643 7626 9050/
!
      call random_number(u)
      do while (u.le.0)                 ! Comment out this block 
        call random_number(u)           ! if your RNG can never
      enddo                             ! return exact zero
      g = c
      u = u+u
      do while (u.lt.1.0)
         g = g + alog2
         u = u+u
      enddo
      u = u-1.0
      if (u.le.p) then
        rexps = g + aa/(bb-u)
        return
      endif
      do
        call random_number(u)
        y = a/(b-u)
        call random_number(up)
        if ((up*hh+dd)*(b-u)**2 .le. exp(-(y+c))) then
          rexps = g+y
          return
        endif
      enddo
      end
SHAR_EOF
fi # end of overwriting check
cd ..
#       End of shell archive
exit 0