dlbeta.f 1.9 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162
  1. *DECK DLBETA
  2. DOUBLE PRECISION FUNCTION DLBETA (A, B)
  3. C***BEGIN PROLOGUE DLBETA
  4. C***PURPOSE Compute the natural logarithm of the complete Beta
  5. C function.
  6. C***LIBRARY SLATEC (FNLIB)
  7. C***CATEGORY C7B
  8. C***TYPE DOUBLE PRECISION (ALBETA-S, DLBETA-D, CLBETA-C)
  9. C***KEYWORDS FNLIB, LOGARITHM OF THE COMPLETE BETA FUNCTION,
  10. C SPECIAL FUNCTIONS
  11. C***AUTHOR Fullerton, W., (LANL)
  12. C***DESCRIPTION
  13. C
  14. C DLBETA(A,B) calculates the double precision natural logarithm of
  15. C the complete beta function for double precision arguments
  16. C A and B.
  17. C
  18. C***REFERENCES (NONE)
  19. C***ROUTINES CALLED D9LGMC, DGAMMA, DLNGAM, DLNREL, XERMSG
  20. C***REVISION HISTORY (YYMMDD)
  21. C 770701 DATE WRITTEN
  22. C 890531 Changed all specific intrinsics to generic. (WRB)
  23. C 890531 REVISION DATE from Version 3.2
  24. C 891214 Prologue converted to Version 4.0 format. (BAB)
  25. C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
  26. C 900727 Added EXTERNAL statement. (WRB)
  27. C***END PROLOGUE DLBETA
  28. DOUBLE PRECISION A, B, P, Q, CORR, SQ2PIL, D9LGMC, DGAMMA, DLNGAM,
  29. 1 DLNREL
  30. EXTERNAL DGAMMA
  31. SAVE SQ2PIL
  32. DATA SQ2PIL / 0.9189385332 0467274178 0329736405 62 D0 /
  33. C***FIRST EXECUTABLE STATEMENT DLBETA
  34. P = MIN (A, B)
  35. Q = MAX (A, B)
  36. C
  37. IF (P .LE. 0.D0) CALL XERMSG ('SLATEC', 'DLBETA',
  38. + 'BOTH ARGUMENTS MUST BE GT ZERO', 1, 2)
  39. C
  40. IF (P.GE.10.D0) GO TO 30
  41. IF (Q.GE.10.D0) GO TO 20
  42. C
  43. C P AND Q ARE SMALL.
  44. C
  45. DLBETA = LOG (DGAMMA(P) * (DGAMMA(Q)/DGAMMA(P+Q)) )
  46. RETURN
  47. C
  48. C P IS SMALL, BUT Q IS BIG.
  49. C
  50. 20 CORR = D9LGMC(Q) - D9LGMC(P+Q)
  51. DLBETA = DLNGAM(P) + CORR + P - P*LOG(P+Q)
  52. 1 + (Q-0.5D0)*DLNREL(-P/(P+Q))
  53. RETURN
  54. C
  55. C P AND Q ARE BIG.
  56. C
  57. 30 CORR = D9LGMC(P) + D9LGMC(Q) - D9LGMC(P+Q)
  58. DLBETA = -0.5D0*LOG(Q) + SQ2PIL + CORR + (P-0.5D0)*LOG(P/(P+Q))
  59. 1 + Q*DLNREL(-P/(P+Q))
  60. RETURN
  61. C
  62. END