dlbeta.f

Go to the documentation of this file.

*DECK DLBETA
      DOUBLE PRECISION FUNCTION dlbeta (A, B)
C***BEGIN PROLOGUE  DLBETA
C***PURPOSE  Compute the natural logarithm of the complete Beta
C            function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7B
C***TYPE      DOUBLE PRECISION (ALBETA-S, DLBETA-D, CLBETA-C)
C***KEYWORDS  FNLIB, LOGARITHM OF THE COMPLETE BETA FUNCTION,
C             SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C DLBETA(A,B) calculates the double precision natural logarithm of
C the complete beta function for double precision arguments
C A and B.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D9LGMC, DGAMMA, DLNGAM, DLNREL, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770701  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   900727  Added EXTERNAL statement.  (WRB)
C***END PROLOGUE  DLBETA
      DOUBLE PRECISION a, b, p, q, corr, sq2pil, d9lgmc, dgamma, dlngam,
     1  dlnrel
      EXTERNAL dgamma
      SAVE sq2pil
      DATA sq2pil / 0.9189385332 0467274178 0329736405 62 d0 /
C***FIRST EXECUTABLE STATEMENT  DLBETA
      p = min(a, b)
      q = max(a, b)
C
      IF (p .LE. 0.d0) CALL xermsg('SLATEC', 'DLBETA',
     +   'BOTH ARGUMENTS MUST BE GT ZERO', 1, 2)
C
      IF (p.GE.10.d0) go to 30
      IF (q.GE.10.d0) go to 20
C
C P AND Q ARE SMALL.
C
      dlbeta = log(dgamma(p) * (dgamma(q)/dgamma(p+q)) )
      RETURN
C
C P IS SMALL, BUT Q IS BIG.
C
 20   corr = d9lgmc(q) - d9lgmc(p+q)
      dlbeta = dlngam(p) + corr + p - p*log(p+q)
     1  + (q-0.5d0)*dlnrel(-p/(p+q))
      RETURN
C
C P AND Q ARE BIG.
C
 30   corr = d9lgmc(p) + d9lgmc(q) - d9lgmc(p+q)
      dlbeta = -0.5d0*log(q) + sq2pil + corr + (p-0.5d0)*log(p/(p+q))
     1  + q*dlnrel(-p/(p+q))
      RETURN
C
      END