      REAL FUNCTION sexpo()
!**********************************************************************C
!                                                                      C
!                                                                      C
!     (STANDARD-)  E X P O N E N T I A L   DISTRIBUTION                C
!                                                                      C
!                                                                      C
!**********************************************************************C
!**********************************************************************C
!                                                                      C
!     FOR DETAILS SEE:                                                 C
!                                                                      C
!               AHRENS, J.H. AND DIETER, U.                            C
!               COMPUTER METHODS FOR SAMPLING FROM THE                 C
!               EXPONENTIAL AND NORMAL DISTRIBUTIONS.                  C
!               COMM. ACM, 15,10 (OCT. 1972), 873 - 882.               C
!                                                                      C
!     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM       C
!     'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)       C
!                                                                      C
!     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   C
!     SUNIF.  The argument IR thus goes away.                          C
!                                                                      C
!**********************************************************************C
!
      DIMENSION q(8)
      EQUIVALENCE (q(1),q1)
!
!     Q(N) = SUM(ALOG(2.0)**K/K!)    K=1,..,N ,      THE HIGHEST N
!     (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
!
      DATA q/.6931472,.9333737,.9888778,.9984959,.9998293,.9999833,     &
     &     .9999986,.9999999/
!
   10 a = 0.0
      u = ranf()
      GO TO 30

   20 a = a + q1
   30 u = u + u
      IF (u.LE.1.0) GO TO 20
   40 u = u - 1.0
      IF (u.GT.q1) GO TO 60
   50 sexpo = a + u
      RETURN

   60 i = 1
      ustar = ranf()
      umin = ustar
   70 ustar = ranf()
      IF (ustar.LT.umin) umin = ustar
   80 i = i + 1
      IF (u.GT.q(i)) GO TO 70
   90 sexpo = a + umin*q1
      RETURN

      END
