double precision function dsinh (x)
c may 1980 edition.  w. fullerton, c3, los alamos scientific lab.
      double precision x, sinhcs(13), sqeps, y, ymax,
     1  d1mach, dcsevl, dexp, dsqrt
      external d1mach, dcsevl, dexp, dsqrt, initds
c
c series for sinh       on the interval  0.          to  1.00000e+00
c                                        with weighted error   7.76e-33
c                                         log weighted error  32.11
c                               significant figures required  31.20
c                                    decimal places required  32.67
c
      data sinhcs(  1) / +.1730421940 4717963167 5883846985 01 d+0     /
      data sinhcs(  2) / +.8759422192 2760477154 9002634544 40 d-1     /
      data sinhcs(  3) / +.1079477774 5671327502 4272706515 79 d-2     /
      data sinhcs(  4) / +.6374849260 7547504815 6855545718 50 d-5     /
      data sinhcs(  5) / +.2202366404 9230530159 1904960195 02 d-7     /
      data sinhcs(  6) / +.4987940180 4158493149 4258072036 61 d-10    /
      data sinhcs(  7) / +.7973053554 1157304814 4114804411 86 d-13    /
      data sinhcs(  8) / +.9473158713 0725443342 9273172266 66 d-16    /
      data sinhcs(  9) / +.8693492050 4480078871 0230986666 66 d-19    /
      data sinhcs( 10) / +.6346939440 3318040457 3973333333 33 d-22    /
      data sinhcs( 11) / +.3774033787 0858485738 6666666666 66 d-25    /
      data sinhcs( 12) / +.1863021371 9570056533 3333333333 33 d-28    /
      data sinhcs( 13) / +.7756843716 6506666666 6666666666 66 d-32    /
c
      data nterms, sqeps, ymax / 0, 2*0.0d0 /
c
      if (nterms.ne.0) go to 10
      nterms = initds (sinhcs, 13, 0.1*sngl(d1mach(3)) )
      sqeps = dsqrt (6.0d0*d1mach(3))
      ymax = 1.0d0/dsqrt(d1mach(3))
c
 10   y = dabs(x)
      if (y.gt.1.d0) go to 20
c
      dsinh = x
      if (y.gt.sqeps) dsinh = x*(1.d0 + dcsevl (2.d0*x*x-1.d0, sinhcs,
     1  nterms) )
      return
c
 20   y = dexp (y)
      if (y.ge.ymax) dsinh = dsign (0.5d0*y, x)
      if (y.lt.ymax) dsinh = dsign (0.5d0*(y-1.0d0/y), x)
c
      return
      end