src/numerical/slatec/fortran/d9lgmc.f

   1 *DECK D9LGMC
   2       DOUBLE PRECISION FUNCTION D9LGMC (X)
   3 C***BEGIN PROLOGUE  D9LGMC
   4 C***SUBSIDIARY
   5 C***PURPOSE  Compute the log Gamma correction factor so that
   6 C            LOG(DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-5.)*LOG(X) - X
   7 C            + D9LGMC(X).
   8 C***LIBRARY   SLATEC (FNLIB)
   9 C***CATEGORY  C7E
  10 C***TYPE      DOUBLE PRECISION (R9LGMC-S, D9LGMC-D, C9LGMC-C)
  11 C***KEYWORDS  COMPLETE GAMMA FUNCTION, CORRECTION TERM, FNLIB,
  12 C             LOG GAMMA, LOGARITHM, SPECIAL FUNCTIONS
  13 C***AUTHOR  Fullerton, W., (LANL)
  14 C***DESCRIPTION
  15 C
  16 C Compute the log gamma correction factor for X .GE. 10. so that
  17 C LOG (DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-.5)*LOG(X) - X + D9lGMC(X)
  18 C
  19 C Series for ALGM       on the interval  0.          to  1.00000E-02
  20 C                                        with weighted error   1.28E-31
  21 C                                         log weighted error  30.89
  22 C                               significant figures required  29.81
  23 C                                    decimal places required  31.48
  24 C
  25 C***REFERENCES  (NONE)
  26 C***ROUTINES CALLED  D1MACH, DCSEVL, INITDS, XERMSG
  27 C***REVISION HISTORY  (YYMMDD)
  28 C   770601  DATE WRITTEN
  29 C   890531  Changed all specific intrinsics to generic.  (WRB)
  30 C   890531  REVISION DATE from Version 3.2
  31 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  32 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
  33 C   900720  Routine changed from user-callable to subsidiary.  (WRB)
  34 C***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./
  54 C***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.
  61 C
  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
  65 C
  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
  70 C
  71  20   D9LGMC = 0.D0
  72       CALL XERMSG ('SLATEC', 'D9LGMC', 'X SO BIG D9LGMC UNDERFLOWS', 2,
  73      +   1)
  74       RETURN
  75 C
  76       END