double precision function dcin (x)
c december 1980 edition.  w. fullerton, bell labs.
      double precision x, cincs(18), eul, xmin, f, g, sinx, absx,
     1  dcsevl, d1mach, dcos, dlog, dsin, dsqrt
      external d1mach, dcos, dcsevl, dlog, dsin, dsqrt, initds
c
c series for cin  on the interval  0.00000e+00 to  1.60000e+01
c                                        with weighted error   2.33e-33
c                                         log weighted error  32.63
c                               significant figures required  31.92
c                                    decimal places required  33.26
c
      data cin cs(  1) /  0.3707450175 0909688741 6548012285 64992d0/
      data cin cs(  2) / -0.0589357489 6364446831 9568643973 63697d0/
      data cin cs(  3) /  0.0053818964 2113569124 0487453262 03340d0/
      data cin cs(  4) / -0.0002986005 2841962135 3195949065 63410d0/
      data cin cs(  5) /  0.0000109557 2575321620 0770310544 67306d0/
      data cin cs(  6) / -0.0000002840 5454877346 6304917271 87731d0/
      data cin cs(  7) /  0.0000000054 6973994875 3849124578 61806d0/
      data cin cs(  8) / -0.0000000000 8124187461 3181570832 77452d0/
      data cin cs(  9) /  0.0000000000 0095868593 1177066090 13181d0/
      data cin cs( 10) / -0.0000000000 0000920266 0043923510 31377d0/
      data cin cs( 11) /  0.0000000000 0000007325 8879990178 95024d0/
      data cin cs( 12) / -0.0000000000 0000000049 1437266758 42909d0/
      data cin cs( 13) /  0.0000000000 0000000000 2815777467 53902d0/
      data cin cs( 14) / -0.0000000000 0000000000 0013939867 88501d0/
      data cin cs( 15) /  0.0000000000 0000000000 0000060224 85646d0/
      data cin cs( 16) / -0.0000000000 0000000000 0000000229 04717d0/
      data cin cs( 17) /  0.0000000000 0000000000 0000000000 77273d0/
      data cin cs( 18) / -0.0000000000 0000000000 0000000000 00233d0/
c
      data eul / 0.5772156649 0153286060 6512090082 40 d0 /
      data ncin, xmin /0, 0.d0/
c
      if (ncin.ne.0) go to 10
      ncin = initds (cincs, 18, 0.1*sngl(d1mach(3)))
      xmin = dsqrt (d1mach(1))
c
 10   dcin = 0.d0
      absx = dabs(x)
      if (absx.ne.0.d0 .and. absx.le.xmin) call seteru (
     1  38hdcin    x so small that cin underflows, 38, 1, 0)
      if (absx.le.xmin) return
c
      if (absx.gt.4.0d0) go to 20
      dcin = dcsevl((x*x-8.d0)*.125d0, cincs, ncin)*x*x
      return
c
 20   call d9sifg (absx, f, g)
      sinx = dsin (absx)
      call erroff
      dcin = -f*sinx + g*dcos(absx) + dlog(absx) + eul
c
      return
      end