home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / slatec-fn / d9lgmc.f next >
Text File  |  1996-09-28  |  3KB  |  77 lines
  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
  77.