function shi (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
      dimension shics(7)
      external csevl, e1, ei, inits, r1mach, sqrt
c
c series for shi  on the interval  0.00000e+00 to  1.40625e-01
c                                        with weighted error   4.67e-20
c                                         log weighted error  19.33
c                               significant figures required  17.07
c                                    decimal places required  19.75
c
      data shi cs(  1) /     0.0078372685 6889009506 95e0/
      data shi cs(  2) /     0.0039227664 9342345639 73e0/
      data shi cs(  3) /     0.0000041346 7878876172 67e0/
      data shi cs(  4) /     0.0000000024 7074803728 83e0/
      data shi cs(  5) /     0.0000000000 0093792955 91e0/
      data shi cs(  6) /     0.0000000000 0000024518 17e0/
      data shi cs(  7) /     0.0000000000 0000000004 67e0/
c
      data nshi, xsml /0, 0.0/
c
      if (nshi.ne.0) go to 10
      nshi = inits (shics, 7, 0.1*r1mach(3))
      xsml = sqrt (r1mach(3))
c
 10   absx = abs(x)
      shi = x
      if (absx.gt.xsml .and. absx.le.0.375) shi = x*(1.0 +
     1  csevl (128.*x*x/9.-1., shics, nshi))
      if (absx.gt.0.375) shi = 0.5*(ei(x) + e1(x))
c
      return
      end