double precision function d9lgmc(x)
c***begin prologue  d9lgmc
c***date written   770601   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  c7e
c***keywords  complete gamma function,correction factor,
c             double precision,gamma function,logarithm,
c             special function
c***author  fullerton, w., (lanl)
c***purpose  computes the  d.p. log gamma correction factor for
c            x .ge. 10. so that dlog(dgamma(x)) = dlog(dsqrt(2*pi)) +
c            (x-5.)*dlog(x) - x + d9lgmc(x)
c***description
c
c compute the log gamma correction factor for x .ge. 10. so that
c dlog (dgamma(x)) = dlog(dsqrt(2*pi)) + (x-.5)*dlog(x) - x + d9lgmc(x)
c
c series for algm       on the interval  0.          to  1.00000e-02
c                                        with weighted error   1.28e-31
c                                         log weighted error  30.89
c                               significant figures required  29.81
c                                    decimal places required  31.48
c***references  (none)
c***routines called  d1mach,dcsevl,initds,xerror
c***end prologue  d9lgmc
      double precision x, algmcs(15), xbig, xmax, dcsevl, d1mach
      data algmcs(  1) / +.1666389480 4518632472 0572965082 2 d+0      /
      data algmcs(  2) / -.1384948176 0675638407 3298605913 5 d-4      /
      data algmcs(  3) / +.9810825646 9247294261 5717154748 7 d-8      /
      data algmcs(  4) / -.1809129475 5724941942 6330626671 9 d-10     /
      data algmcs(  5) / +.6221098041 8926052271 2601554341 6 d-13     /
      data algmcs(  6) / -.3399615005 4177219443 0333059966 6 d-15     /
      data algmcs(  7) / +.2683181998 4826987489 5753884666 6 d-17     /
      data algmcs(  8) / -.2868042435 3346432841 4462239999 9 d-19     /
      data algmcs(  9) / +.3962837061 0464348036 7930666666 6 d-21     /
      data algmcs( 10) / -.6831888753 9857668701 1199999999 9 d-23     /
      data algmcs( 11) / +.1429227355 9424981475 7333333333 3 d-24     /
      data algmcs( 12) / -.3547598158 1010705471 9999999999 9 d-26     /
      data algmcs( 13) / +.1025680058 0104709120 0000000000 0 d-27     /
      data algmcs( 14) / -.3401102254 3167487999 9999999999 9 d-29     /
      data algmcs( 15) / +.1276642195 6300629333 3333333333 3 d-30     /
      data nalgm, xbig, xmax / 0, 2*0.d0 /
c***first executable statement  d9lgmc
      if (nalgm.ne.0) go to 10
      nalgm = initds (algmcs, 15, sngl(d1mach(3)) )
      xbig = 1.0d0/dsqrt(d1mach(3))
      xmax = dexp (dmin1(dlog(d1mach(2)/12.d0), -dlog(12.d0*d1mach(1))))
c
 10   if (x.lt.10.d0) call xerror ( 'd9lgmc  x must be ge 10', 23, 1, 2)
      if (x.ge.xmax) go to 20
c
      d9lgmc = 1.d0/(12.d0*x)
      if (x.lt.xbig) d9lgmc = dcsevl (2.0d0*(10.d0/x)**2-1.d0, algmcs,
     1  nalgm) / x
      return
c
 20   d9lgmc = 0.d0
      call xerror ( 'd9lgmc  x so big d9lgmc underflows', 34, 2, 1)
      return
c
      end