GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
albeta.f
Go to the documentation of this file.
1 *DECK ALBETA
2  FUNCTION albeta (A, B)
3 C***BEGIN PROLOGUE ALBETA
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 SINGLE 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 ALBETA computes the natural log of the complete beta function.
15 C
16 C Input Parameters:
17 C A real and positive
18 C B real and positive
19 C
20 C***REFERENCES (NONE)
21 C***ROUTINES CALLED ALNGAM, ALNREL, GAMMA, R9LGMC, XERMSG
22 C***REVISION HISTORY (YYMMDD)
23 C 770701 DATE WRITTEN
24 C 890531 Changed all specific intrinsics to generic. (WRB)
25 C 890531 REVISION DATE from Version 3.2
26 C 891214 Prologue converted to Version 4.0 format. (BAB)
27 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
28 C 900326 Removed duplicate information from DESCRIPTION section.
29 C (WRB)
30 C 900727 Added EXTERNAL statement. (WRB)
31 C***END PROLOGUE ALBETA
32  EXTERNAL gamma
33  SAVE sq2pil
34  DATA sq2pil / 0.9189385332 0467274 e0 /
35 C***FIRST EXECUTABLE STATEMENT ALBETA
36  p = min(a, b)
37  q = max(a, b)
38 C
39  IF (p .LE. 0.0) CALL xermsg('SLATEC', 'ALBETA',
40  + 'BOTH ARGUMENTS MUST BE GT ZERO', 1, 2)
41  IF (p.GE.10.0) go to 30
42  IF (q.GE.10.0) go to 20
43 C
44 C P AND Q ARE SMALL.
45 C
46  albeta = log(gamma(p) * (gamma(q)/gamma(p+q)) )
47  RETURN
48 C
49 C P IS SMALL, BUT Q IS BIG.
50 C
51  20 corr = r9lgmc(q) - r9lgmc(p+q)
52  albeta = alngam(p) + corr + p - p*log(p+q) +
53  1 (q-0.5)*alnrel(-p/(p+q))
54  RETURN
55 C
56 C P AND Q ARE BIG.
57 C
58  30 corr = r9lgmc(p) + r9lgmc(q) - r9lgmc(p+q)
59  albeta = -0.5*log(q) + sq2pil + corr + (p-0.5)*log(p/(p+q))
60  1 + q*alnrel(-p/(p+q))
61  RETURN
62 C
63  END