double precision function dexprl (x)
  c july 1980 edition.  w. fullerton, c3, los alamos scientific lab.
  c
  c evaluate  exprel(x) = (exp(x) - 1.0) / x.   for small abs(x) the
  c taylor series is used.  if x is negative the reflection formula
  c         exprel(x) = exp(x) * exprel(abs(x))
  c may be used  this reflection formula will be of use when the
  c evaluation for small abs(x) is done by chebyshev series rather than
  c taylor series.
  c
        double precision x, absx, alneps, xbnd, xln, xn,  d1mach,
       1  dexp, dlog
        external d1mach
        data nterms, xbnd / 0, 0.0d0 /
  c
        if (nterms.ne.0) go to 10
        alneps = dlog(d1mach(3))
        xn = 3.72d0 - 0.3d0*alneps
        xln = dlog((xn+1.0d0)/1.36d0)
        nterms = xn - (xn*xln+alneps)/(xln+1.36d0) + 1.5d0
        xbnd = d1mach(3)
  c
   10   absx = dabs(x)
        if (absx.le.0.5d0) go to 20
        call entsrc (irold, 1)
        dexprl = (dexp(x)-1.0d0)/x
        call entsrc (ir, irold)
        return
  c
   20   dexprl = 1.0d0
        if (absx.lt.xbnd) return
  c
        dexprl = 0.0d0
        do 30 i=1,nterms
          dexprl = 1.0d0 + dexprl*x/dble(float(nterms+2-i))
   30   continue
  c
        return
        end