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 / dgamma.f < prev    next >
Text File  |  1996-09-28  |  7KB  |  154 lines
  1. *DECK DGAMMA
  2.       DOUBLE PRECISION FUNCTION DGAMMA (X)
  3. C***BEGIN PROLOGUE  DGAMMA
  4. C***PURPOSE  Compute the complete Gamma function.
  5. C***LIBRARY   SLATEC (FNLIB)
  6. C***CATEGORY  C7A
  7. C***TYPE      DOUBLE PRECISION (GAMMA-S, DGAMMA-D, CGAMMA-C)
  8. C***KEYWORDS  COMPLETE GAMMA FUNCTION, FNLIB, SPECIAL FUNCTIONS
  9. C***AUTHOR  Fullerton, W., (LANL)
  10. C***DESCRIPTION
  11. C
  12. C DGAMMA(X) calculates the double precision complete Gamma function
  13. C for double precision argument X.
  14. C
  15. C Series for GAM        on the interval  0.          to  1.00000E+00
  16. C                                        with weighted error   5.79E-32
  17. C                                         log weighted error  31.24
  18. C                               significant figures required  30.00
  19. C                                    decimal places required  32.05
  20. C
  21. C***REFERENCES  (NONE)
  22. C***ROUTINES CALLED  D1MACH, D9LGMC, DCSEVL, DGAMLM, INITDS, XERMSG
  23. C***REVISION HISTORY  (YYMMDD)
  24. C   770601  DATE WRITTEN
  25. C   890531  Changed all specific intrinsics to generic.  (WRB)
  26. C   890911  Removed unnecessary intrinsics.  (WRB)
  27. C   890911  REVISION DATE from Version 3.2
  28. C   891214  Prologue converted to Version 4.0 format.  (BAB)
  29. C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
  30. C   920618  Removed space from variable name.  (RWC, WRB)
  31. C***END PROLOGUE  DGAMMA
  32.       DOUBLE PRECISION X, GAMCS(42), DXREL, PI, SINPIY, SQ2PIL, XMAX,
  33.      1  XMIN, Y, D9LGMC, DCSEVL, D1MACH
  34.       LOGICAL FIRST
  35. C
  36.       SAVE GAMCS, PI, SQ2PIL, NGAM, XMIN, XMAX, DXREL, FIRST
  37.       DATA GAMCS(  1) / +.8571195590 9893314219 2006239994 2 D-2      /
  38.       DATA GAMCS(  2) / +.4415381324 8410067571 9131577165 2 D-2      /
  39.       DATA GAMCS(  3) / +.5685043681 5993633786 3266458878 9 D-1      /
  40.       DATA GAMCS(  4) / -.4219835396 4185605010 1250018662 4 D-2      /
  41.       DATA GAMCS(  5) / +.1326808181 2124602205 8400679635 2 D-2      /
  42.       DATA GAMCS(  6) / -.1893024529 7988804325 2394702388 6 D-3      /
  43.       DATA GAMCS(  7) / +.3606925327 4412452565 7808221722 5 D-4      /
  44.       DATA GAMCS(  8) / -.6056761904 4608642184 8554829036 5 D-5      /
  45.       DATA GAMCS(  9) / +.1055829546 3022833447 3182350909 3 D-5      /
  46.       DATA GAMCS( 10) / -.1811967365 5423840482 9185589116 6 D-6      /
  47.       DATA GAMCS( 11) / +.3117724964 7153222777 9025459316 9 D-7      /
  48.       DATA GAMCS( 12) / -.5354219639 0196871408 7408102434 7 D-8      /
  49.       DATA GAMCS( 13) / +.9193275519 8595889468 8778682594 0 D-9      /
  50.       DATA GAMCS( 14) / -.1577941280 2883397617 6742327395 3 D-9      /
  51.       DATA GAMCS( 15) / +.2707980622 9349545432 6654043308 9 D-10     /
  52.       DATA GAMCS( 16) / -.4646818653 8257301440 8166105893 3 D-11     /
  53.       DATA GAMCS( 17) / +.7973350192 0074196564 6076717535 9 D-12     /
  54.       DATA GAMCS( 18) / -.1368078209 8309160257 9949917230 9 D-12     /
  55.       DATA GAMCS( 19) / +.2347319486 5638006572 3347177168 8 D-13     /
  56.       DATA GAMCS( 20) / -.4027432614 9490669327 6657053469 9 D-14     /
  57.       DATA GAMCS( 21) / +.6910051747 3721009121 3833697525 7 D-15     /
  58.       DATA GAMCS( 22) / -.1185584500 2219929070 5238712619 2 D-15     /
  59.       DATA GAMCS( 23) / +.2034148542 4963739552 0102605193 2 D-16     /
  60.       DATA GAMCS( 24) / -.3490054341 7174058492 7401294910 8 D-17     /
  61.       DATA GAMCS( 25) / +.5987993856 4853055671 3505106602 6 D-18     /
  62.       DATA GAMCS( 26) / -.1027378057 8722280744 9006977843 1 D-18     /
  63.       DATA GAMCS( 27) / +.1762702816 0605298249 4275966074 8 D-19     /
  64.       DATA GAMCS( 28) / -.3024320653 7353062609 5877211204 2 D-20     /
  65.       DATA GAMCS( 29) / +.5188914660 2183978397 1783355050 6 D-21     /
  66.       DATA GAMCS( 30) / -.8902770842 4565766924 4925160106 6 D-22     /
  67.       DATA GAMCS( 31) / +.1527474068 4933426022 7459689130 6 D-22     /
  68.       DATA GAMCS( 32) / -.2620731256 1873629002 5732833279 9 D-23     /
  69.       DATA GAMCS( 33) / +.4496464047 8305386703 3104657066 6 D-24     /
  70.       DATA GAMCS( 34) / -.7714712731 3368779117 0390152533 3 D-25     /
  71.       DATA GAMCS( 35) / +.1323635453 1260440364 8657271466 6 D-25     /
  72.       DATA GAMCS( 36) / -.2270999412 9429288167 0231381333 3 D-26     /
  73.       DATA GAMCS( 37) / +.3896418998 0039914493 2081663999 9 D-27     /
  74.       DATA GAMCS( 38) / -.6685198115 1259533277 9212799999 9 D-28     /
  75.       DATA GAMCS( 39) / +.1146998663 1400243843 4761386666 6 D-28     /
  76.       DATA GAMCS( 40) / -.1967938586 3451346772 9510399999 9 D-29     /
  77.       DATA GAMCS( 41) / +.3376448816 5853380903 3489066666 6 D-30     /
  78.       DATA GAMCS( 42) / -.5793070335 7821357846 2549333333 3 D-31     /
  79.       DATA PI / 3.1415926535 8979323846 2643383279 50 D0 /
  80.       DATA SQ2PIL / 0.9189385332 0467274178 0329736405 62 D0 /
  81.       DATA FIRST /.TRUE./
  82. C***FIRST EXECUTABLE STATEMENT  DGAMMA
  83.       IF (FIRST) THEN
  84.          NGAM = INITDS (GAMCS, 42, 0.1*REAL(D1MACH(3)) )
  85. C
  86.          CALL DGAMLM (XMIN, XMAX)
  87.          DXREL = SQRT(D1MACH(4))
  88.       ENDIF
  89.       FIRST = .FALSE.
  90. C
  91.       Y = ABS(X)
  92.       IF (Y.GT.10.D0) GO TO 50
  93. C
  94. C COMPUTE GAMMA(X) FOR -XBND .LE. X .LE. XBND.  REDUCE INTERVAL AND FIND
  95. C GAMMA(1+Y) FOR 0.0 .LE. Y .LT. 1.0 FIRST OF ALL.
  96. C
  97.       N = X
  98.       IF (X.LT.0.D0) N = N - 1
  99.       Y = X - N
  100.       N = N - 1
  101.       DGAMMA = 0.9375D0 + DCSEVL (2.D0*Y-1.D0, GAMCS, NGAM)
  102.       IF (N.EQ.0) RETURN
  103. C
  104.       IF (N.GT.0) GO TO 30
  105. C
  106. C COMPUTE GAMMA(X) FOR X .LT. 1.0
  107. C
  108.       N = -N
  109.       IF (X .EQ. 0.D0) CALL XERMSG ('SLATEC', 'DGAMMA', 'X IS 0', 4, 2)
  110.       IF (X .LT. 0.0 .AND. X+N-2 .EQ. 0.D0) CALL XERMSG ('SLATEC',
  111.      +   'DGAMMA', 'X IS A NEGATIVE INTEGER', 4, 2)
  112.       IF (X .LT. (-0.5D0) .AND. ABS((X-AINT(X-0.5D0))/X) .LT. DXREL)
  113.      +   CALL XERMSG ('SLATEC', 'DGAMMA',
  114.      +   'ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER',
  115.      +   1, 1)
  116. C
  117.       DO 20 I=1,N
  118.         DGAMMA = DGAMMA/(X+I-1 )
  119.  20   CONTINUE
  120.       RETURN
  121. C
  122. C GAMMA(X) FOR X .GE. 2.0 AND X .LE. 10.0
  123. C
  124.  30   DO 40 I=1,N
  125.         DGAMMA = (Y+I) * DGAMMA
  126.  40   CONTINUE
  127.       RETURN
  128. C
  129. C GAMMA(X) FOR ABS(X) .GT. 10.0.  RECALL Y = ABS(X).
  130. C
  131.  50   IF (X .GT. XMAX) CALL XERMSG ('SLATEC', 'DGAMMA',
  132.      +   'X SO BIG GAMMA OVERFLOWS', 3, 2)
  133. C
  134.       DGAMMA = 0.D0
  135.       IF (X .LT. XMIN) CALL XERMSG ('SLATEC', 'DGAMMA',
  136.      +   'X SO SMALL GAMMA UNDERFLOWS', 2, 1)
  137.       IF (X.LT.XMIN) RETURN
  138. C
  139.       DGAMMA = EXP ((Y-0.5D0)*LOG(Y) - Y + SQ2PIL + D9LGMC(Y) )
  140.       IF (X.GT.0.D0) RETURN
  141. C
  142.       IF (ABS((X-AINT(X-0.5D0))/X) .LT. DXREL) CALL XERMSG ('SLATEC',
  143.      +   'DGAMMA',
  144.      +   'ANSWER LT HALF PRECISION, X TOO NEAR NEGATIVE INTEGER', 1, 1)
  145. C
  146.       SINPIY = SIN (PI*Y)
  147.       IF (SINPIY .EQ. 0.D0) CALL XERMSG ('SLATEC', 'DGAMMA',
  148.      +   'X IS A NEGATIVE INTEGER', 4, 2)
  149. C
  150.       DGAMMA = -PI/(Y*SINPIY*DGAMMA)
  151. C
  152.       RETURN
  153.       END
  154.