double precision function dshi (x)
c december 1980 edition.  w. fullerton, bell labs.
c
c evaluate the hyperbolic sine integral
c            shi = integral from 0 to x of  sinh(t)/t dt.
c
      double precision x, shics(10), xsml, absx,
     1  d1mach, dcsevl, de1, dei, dsqrt
c
c series for shi  on the interval  0.00000e+00 to  1.40625e-01
c                                        with weighted error   7.11e-32
c                                         log weighted error  31.15
c                               significant figures required  28.89
c                                    decimal places required  31.65
c
      data shi cs(  1) /     0.0078372685 6889009506 9520098431 7332d0/
      data shi cs(  2) /     0.0039227664 9342345639 7269757442 7225d0/
      data shi cs(  3) /     0.0000041346 7878876172 6674674790 8275d0/
      data shi cs(  4) /     0.0000000024 7074803728 8274213514 5302d0/
      data shi cs(  5) /     0.0000000000 0093792955 9076363045 7157d0/
      data shi cs(  6) /     0.0000000000 0000024518 1701952086 7353d0/
      data shi cs(  7) /     0.0000000000 0000000004 6741615525 7592d0/
      data shi cs(  8) /     0.0000000000 0000000000 0006780307 2389d0/
      data shi cs(  9) /     0.0000000000 0000000000 0000000773 1289d0/
      data shi cs( 10) /     0.0000000000 0000000000 0000000000 0711d0/
c
      data nshi, xsml /0, 0.0d0/
c
      if (nshi.ne.0) go to 10
      nshi = initds (shics, 10, 0.1*sngl(d1mach(3)))
      xsml = dsqrt (d1mach(3))
c
 10   absx = dabs(x)
c
      dshi = x
      if (absx.gt.xsml .and. absx.le.0.375d0) dshi = x*(1.d0 +
     1  dcsevl (128.d0*x*x/9.d0-1.d0, shics, nshi))
      if (absx.gt.0.375d0) dshi = 0.5d0*(dei(x) + de1(x))
c
      return
      end