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