d9lgmc.f

Go to the documentation of this file.

*DECK D9LGMC
      DOUBLE PRECISION FUNCTION d9lgmc (X)
C***BEGIN PROLOGUE  D9LGMC
C***SUBSIDIARY
C***PURPOSE  Compute the log Gamma correction factor so that
C            LOG(DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-5.)*LOG(X) - X
C            + D9LGMC(X).
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7E
C***TYPE      DOUBLE PRECISION (R9LGMC-S, D9LGMC-D, C9LGMC-C)
C***KEYWORDS  COMPLETE GAMMA FUNCTION, CORRECTION TERM, FNLIB,
C             LOG GAMMA, LOGARITHM, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C Compute the log gamma correction factor for X .GE. 10. so that
C LOG (DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-.5)*LOG(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
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, DCSEVL, INITDS, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900720  Routine changed from user-callable to subsidiary.  (WRB)
C***END PROLOGUE  D9LGMC
      DOUBLE PRECISION x, algmcs(15), xbig, xmax, dcsevl, d1mach
      LOGICAL first
      SAVE algmcs, nalgm, xbig, xmax, first
      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 first /.true./
C***FIRST EXECUTABLE STATEMENT  D9LGMC
      IF (first) THEN
         nalgm = initds(algmcs, 15, real(d1mach(3)) )
         xbig = 1.0d0/sqrt(d1mach(3))
         xmax = exp(min(log(d1mach(2)/12.d0), -log(12.d0*d1mach(1))))
      ENDIF
      first = .false.
C
      IF (x .LT. 10.d0) CALL xermsg ('SLATEC', 'D9LGMC',
     +   'X MUST BE GE 10', 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 xermsg ('SLATEC', 'D9LGMC', 'X SO BIG D9LGMC UNDERFLOWS', 2,
     +   1)
      RETURN
C
      END