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
  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