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
alngam.f
Go to the documentation of this file.
1 *DECK ALNGAM
2  FUNCTION alngam (X)
3 C***BEGIN PROLOGUE ALNGAM
4 C***PURPOSE Compute the logarithm of the absolute value of the Gamma
5 C function.
6 C***LIBRARY SLATEC (FNLIB)
7 C***CATEGORY C7A
8 C***TYPE SINGLE PRECISION (ALNGAM-S, DLNGAM-D, CLNGAM-C)
9 C***KEYWORDS ABSOLUTE VALUE, COMPLETE GAMMA FUNCTION, FNLIB, LOGARITHM,
10 C SPECIAL FUNCTIONS
11 C***AUTHOR Fullerton, W., (LANL)
12 C***DESCRIPTION
13 C
14 C ALNGAM(X) computes the logarithm of the absolute value of the
15 C gamma function at X.
16 C
17 C***REFERENCES (NONE)
18 C***ROUTINES CALLED GAMMA, R1MACH, R9LGMC, XERMSG
19 C***REVISION HISTORY (YYMMDD)
20 C 770601 DATE WRITTEN
21 C 890531 Changed all specific intrinsics to generic. (WRB)
22 C 890531 REVISION DATE from Version 3.2
23 C 891214 Prologue converted to Version 4.0 format. (BAB)
24 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
25 C 900326 Removed duplicate information from DESCRIPTION section.
26 C (WRB)
27 C 900727 Added EXTERNAL statement. (WRB)
28 C***END PROLOGUE ALNGAM
29  LOGICAL first
30  EXTERNAL gamma
31  SAVE sq2pil, sqpi2l, pi, xmax, dxrel, first
32  DATA sq2pil / 0.9189385332 0467274e0/
33  DATA sqpi2l / 0.2257913526 4472743e0/
34  DATA pi / 3.1415926535 8979324e0/
35  DATA first /.true./
36 C***FIRST EXECUTABLE STATEMENT ALNGAM
37  IF (first) THEN
38  xmax = r1mach(2)/log(r1mach(2))
39  dxrel = sqrt(r1mach(4))
40  ENDIF
41  first = .false.
42 C
43  y = abs(x)
44  IF (y.GT.10.0) go to 20
45 C
46 C LOG (ABS (GAMMA(X))) FOR ABS(X) .LE. 10.0
47 C
48  alngam = log(abs(gamma(x)))
49  RETURN
50 C
51 C LOG (ABS (GAMMA(X))) FOR ABS(X) .GT. 10.0
52 C
53  20 IF (y .GT. xmax) CALL xermsg('SLATEC', 'ALNGAM',
54  + 'ABS(X) SO BIG ALNGAM OVERFLOWS', 2, 2)
55 C
56  IF (x.GT.0.) alngam = sq2pil + (x-0.5)*log(x) - x + r9lgmc(y)
57  IF (x.GT.0.) RETURN
58 C
59  sinpiy = abs(sin(pi*y))
60  IF (sinpiy .EQ. 0.) CALL xermsg('SLATEC', 'ALNGAM',
61  + 'X IS A NEGATIVE INTEGER', 3, 2)
62 C
63  IF (abs((x-aint(x-0.5))/x) .LT. dxrel) CALL xermsg('SLATEC',
64  + 'ALNGAM', 'ANSWER LT HALF PRECISION BECAUSE X TOO NEAR ' //
65  + 'NEGATIVE INTEGER', 1, 1)
66 C
67  alngam = sqpi2l + (x-0.5)*log(y) - x - log(sinpiy) - r9lgmc(y)
68  RETURN
69 C
70  END