1*DECK D9LGMC
2      DOUBLE PRECISION FUNCTION D9LGMC (X)
3C***BEGIN PROLOGUE  D9LGMC
4C***SUBSIDIARY
5C***PURPOSE  Compute the log Gamma correction factor so that
6C            LOG(DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-5.)*LOG(X) - X
7C            + D9LGMC(X).
8C***LIBRARY   SLATEC (FNLIB)
9C***CATEGORY  C7E
10C***TYPE      DOUBLE PRECISION (R9LGMC-S, D9LGMC-D, C9LGMC-C)
11C***KEYWORDS  COMPLETE GAMMA FUNCTION, CORRECTION TERM, FNLIB,
12C             LOG GAMMA, LOGARITHM, SPECIAL FUNCTIONS
13C***AUTHOR  Fullerton, W., (LANL)
14C***DESCRIPTION
15C
16C Compute the log gamma correction factor for X .GE. 10. so that
17C LOG (DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-.5)*LOG(X) - X + D9lGMC(X)
18C
19C Series for ALGM       on the interval  0.          to  1.00000E-02
20C                                        with weighted error   1.28E-31
21C                                         log weighted error  30.89
22C                               significant figures required  29.81
23C                                    decimal places required  31.48
24C
25C***REFERENCES  (NONE)
26C***ROUTINES CALLED  D1MACH, DCSEVL, INITDS, XERMSG
27C***REVISION HISTORY  (YYMMDD)
28C   770601  DATE WRITTEN
29C   890531  Changed all specific intrinsics to generic.  (WRB)
30C   890531  REVISION DATE from Version 3.2
31C   891214  Prologue converted to Version 4.0 format.  (BAB)
32C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
33C   900720  Routine changed from user-callable to subsidiary.  (WRB)
34C***END PROLOGUE  D9LGMC
35      DOUBLE PRECISION X, ALGMCS(15), XBIG, XMAX, DCSEVL, D1MACH
36      LOGICAL FIRST
37      SAVE ALGMCS, NALGM, XBIG, XMAX, FIRST
38      DATA ALGMCS(  1) / +.1666389480 4518632472 0572965082 2 D+0      /
39      DATA ALGMCS(  2) / -.1384948176 0675638407 3298605913 5 D-4      /
40      DATA ALGMCS(  3) / +.9810825646 9247294261 5717154748 7 D-8      /
41      DATA ALGMCS(  4) / -.1809129475 5724941942 6330626671 9 D-10     /
42      DATA ALGMCS(  5) / +.6221098041 8926052271 2601554341 6 D-13     /
43      DATA ALGMCS(  6) / -.3399615005 4177219443 0333059966 6 D-15     /
44      DATA ALGMCS(  7) / +.2683181998 4826987489 5753884666 6 D-17     /
45      DATA ALGMCS(  8) / -.2868042435 3346432841 4462239999 9 D-19     /
46      DATA ALGMCS(  9) / +.3962837061 0464348036 7930666666 6 D-21     /
47      DATA ALGMCS( 10) / -.6831888753 9857668701 1199999999 9 D-23     /
48      DATA ALGMCS( 11) / +.1429227355 9424981475 7333333333 3 D-24     /
49      DATA ALGMCS( 12) / -.3547598158 1010705471 9999999999 9 D-26     /
50      DATA ALGMCS( 13) / +.1025680058 0104709120 0000000000 0 D-27     /
51      DATA ALGMCS( 14) / -.3401102254 3167487999 9999999999 9 D-29     /
52      DATA ALGMCS( 15) / +.1276642195 6300629333 3333333333 3 D-30     /
53      DATA FIRST /.TRUE./
54C***FIRST EXECUTABLE STATEMENT  D9LGMC
55      IF (FIRST) THEN
56         NALGM = INITDS (ALGMCS, 15, REAL(D1MACH(3)) )
57         XBIG = 1.0D0/SQRT(D1MACH(3))
58         XMAX = EXP (MIN(LOG(D1MACH(2)/12.D0), -LOG(12.D0*D1MACH(1))))
59      ENDIF
60      FIRST = .FALSE.
61C
62      IF (X .LT. 10.D0) CALL XERMSG ('SLATEC', 'D9LGMC',
63     +   'X MUST BE GE 10', 1, 2)
64      IF (X.GE.XMAX) GO TO 20
65C
66      D9LGMC = 1.D0/(12.D0*X)
67      IF (X.LT.XBIG) D9LGMC = DCSEVL (2.0D0*(10.D0/X)**2-1.D0, ALGMCS,
68     1  NALGM) / X
69      RETURN
70C
71 20   D9LGMC = 0.D0
72      CALL XERMSG ('SLATEC', 'D9LGMC', 'X SO BIG D9LGMC UNDERFLOWS', 2,
73     +   1)
74      RETURN
75C
76      END
77