FUNCTION REXPU() RESULT(fn_val)
!
!     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.
!

IMPLICIT NONE
REAL :: fn_val

REAL, PARAMETER :: alog2 = 0.6931471805599453,  a =  5.7133631526454228,  &
                       b = 3.4142135623730950,  c = -1.6734053240284925,  &
                       p = 0.9802581434685472, aa =  5.6005707569738080,  &
                      bb = 3.3468106480569850, hh =  0.0026106723602095,  &
                      dd = 0.0857864376269050
REAL :: g, u, up, y
!
!       Comment out the following lines if your RNG can never return 0.0
DO
  call random_number(u)
  if (u > 0.0) EXIT               ! Zero-protector
END DO

g = c
DO
  u = u + u
  if (u >= 1.0) EXIT
  g = g + alog2
END DO

u = u - 1.0
if (u <= p) THEN
  fn_val = g + aa/(bb-u)
  return
END IF
DO
  call random_number(u)
  y = a/(b-u)
  call random_number(up)
  if ((up*hh + dd)*(b-u)**2 <= exp(-(y+c))) EXIT
END DO

fn_val = g + y

return
end FUNCTION REXPU


FUNCTION REXPS() RESULT(fn_val)
!
!     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.
!
IMPLICIT NONE
REAL :: fn_val

REAL, PARAMETER :: alog2 = 0.6931471805599453,  a =  5.7133631526454228,  &
                       b = 3.4142135623730950,  c = -1.6734053240284925,  &
                       p = 0.9802581434685472, aa =  5.6005707569738080,  &
                      bb = 3.3468106480569850, hh =  0.0026106723602095,  &
                      dd = 0.0857864376269050
REAL :: g, u, up, y
!
!       Comment out the following lines if your RNG can never return 0.0
DO
  call random_number(u)
  if (u > 0.0) EXIT               ! Zero-protector
END DO

g = c
u = u + u
do while (u < 1.0)
   g = g + alog2
   u = u + u
end do
u = u - 1.0
if (u <= p) then
  fn_val = g + aa/(bb-u)
  return
end if

do
  call random_number(u)
  y = a/(b-u)
  call random_number(up)
  if ((up*hh + dd)*(b-u)**2 <= exp(-(y+c))) then
    fn_val = g + y
    return
  end if
end do

RETURN
end FUNCTION REXPS