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
gamma.f
Go to the documentation of this file.
1 *DECK GAMMA
2  FUNCTION gamma (X)
3 C***BEGIN PROLOGUE GAMMA
4 C***PURPOSE Compute the complete Gamma function.
5 C***LIBRARY SLATEC (FNLIB)
6 C***CATEGORY C7A
7 C***TYPE SINGLE 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 GAMMA computes the gamma function at X, where X is not 0, -1, -2, ....
13 C GAMMA and X are single precision.
14 C
15 C***REFERENCES (NONE)
16 C***ROUTINES CALLED CSEVL, GAMLIM, INITS, R1MACH, R9LGMC, XERMSG
17 C***REVISION HISTORY (YYMMDD)
18 C 770601 DATE WRITTEN
19 C 890531 Changed all specific intrinsics to generic. (WRB)
20 C 890531 REVISION DATE from Version 3.2
21 C 891214 Prologue converted to Version 4.0 format. (BAB)
22 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
23 C***END PROLOGUE GAMMA
24  dimension gcs(23)
25  LOGICAL first
26  SAVE gcs, pi, sq2pil, ngcs, xmin, xmax, dxrel, first
27  DATA gcs( 1) / .0085711955 90989331e0/
28  DATA gcs( 2) / .0044153813 24841007e0/
29  DATA gcs( 3) / .0568504368 1599363e0/
30  DATA gcs( 4) /-.0042198353 96418561e0/
31  DATA gcs( 5) / .0013268081 81212460e0/
32  DATA gcs( 6) /-.0001893024 529798880e0/
33  DATA gcs( 7) / .0000360692 532744124e0/
34  DATA gcs( 8) /-.0000060567 619044608e0/
35  DATA gcs( 9) / .0000010558 295463022e0/
36  DATA gcs(10) /-.0000001811 967365542e0/
37  DATA gcs(11) / .0000000311 772496471e0/
38  DATA gcs(12) /-.0000000053 542196390e0/
39  DATA gcs(13) / .0000000009 193275519e0/
40  DATA gcs(14) /-.0000000001 577941280e0/
41  DATA gcs(15) / .0000000000 270798062e0/
42  DATA gcs(16) /-.0000000000 046468186e0/
43  DATA gcs(17) / .0000000000 007973350e0/
44  DATA gcs(18) /-.0000000000 001368078e0/
45  DATA gcs(19) / .0000000000 000234731e0/
46  DATA gcs(20) /-.0000000000 000040274e0/
47  DATA gcs(21) / .0000000000 000006910e0/
48  DATA gcs(22) /-.0000000000 000001185e0/
49  DATA gcs(23) / .0000000000 000000203e0/
50  DATA pi /3.14159 26535 89793 24e0/
51 C SQ2PIL IS LOG (SQRT (2.*PI) )
52  DATA sq2pil /0.91893 85332 04672 74e0/
53  DATA first /.true./
54 C
55 C LANL DEPENDENT CODE REMOVED 81.02.04
56 C
57 C***FIRST EXECUTABLE STATEMENT GAMMA
58  IF (first) THEN
59 C
60 C ---------------------------------------------------------------------
61 C INITIALIZE. FIND LEGAL BOUNDS FOR X, AND DETERMINE THE NUMBER OF
62 C TERMS IN THE SERIES REQUIRED TO ATTAIN AN ACCURACY TEN TIMES BETTER
63 C THAN MACHINE PRECISION.
64 C
65  ngcs = inits(gcs, 23, 0.1*r1mach(3))
66 C
67  CALL gamlim(xmin, xmax)
68  dxrel = sqrt(r1mach(4))
69 C
70 C ---------------------------------------------------------------------
71 C FINISH INITIALIZATION. START EVALUATING GAMMA(X).
72 C
73  ENDIF
74  first = .false.
75 C
76  y = abs(x)
77  IF (y.GT.10.0) go to 50
78 C
79 C COMPUTE GAMMA(X) FOR ABS(X) .LE. 10.0. REDUCE INTERVAL AND
80 C FIND GAMMA(1+Y) FOR 0. .LE. Y .LT. 1. FIRST OF ALL.
81 C
82  n = x
83  IF (x.LT.0.) n = n - 1
84  y = x - n
85  n = n - 1
86  gamma = 0.9375 + csevl(2.*y-1., gcs, ngcs)
87  IF (n.EQ.0) RETURN
88 C
89  IF (n.GT.0) go to 30
90 C
91 C COMPUTE GAMMA(X) FOR X .LT. 1.
92 C
93  n = -n
94  IF (x .EQ. 0.) CALL xermsg('SLATEC', 'GAMMA', 'X IS 0', 4, 2)
95  IF (x .LT. 0. .AND. x+n-2 .EQ. 0.) CALL xermsg('SLATEC', 'GAMMA'
96  1, 'X IS A NEGATIVE INTEGER', 4, 2)
97  IF (x .LT. (-0.5) .AND. abs((x-aint(x-0.5))/x) .LT. dxrel) CALL
98  1xermsg( 'SLATEC', 'GAMMA',
99  2'ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER'
100  3, 1, 1)
101 C
102  DO 20 i=1,n
103  gamma = gamma / (x+i-1)
104  20 CONTINUE
105  RETURN
106 C
107 C GAMMA(X) FOR X .GE. 2.
108 C
109  30 DO 40 i=1,n
110  gamma = (y+i)*gamma
111  40 CONTINUE
112  RETURN
113 C
114 C COMPUTE GAMMA(X) FOR ABS(X) .GT. 10.0. RECALL Y = ABS(X).
115 C
116  50 IF (x .GT. xmax) CALL xermsg('SLATEC', 'GAMMA',
117  + 'X SO BIG GAMMA OVERFLOWS', 3, 2)
118 C
119  gamma = 0.
120  IF (x .LT. xmin) CALL xermsg('SLATEC', 'GAMMA',
121  + 'X SO SMALL GAMMA UNDERFLOWS', 2, 1)
122  IF (x.LT.xmin) RETURN
123 C
124  gamma = exp((y-0.5)*log(y) - y + sq2pil + r9lgmc(y) )
125  IF (x.GT.0.) RETURN
126 C
127  IF (abs((x-aint(x-0.5))/x) .LT. dxrel) CALL xermsg('SLATEC',
128  + 'GAMMA',
129  + 'ANSWER LT HALF PRECISION, X TOO NEAR NEGATIVE INTEGER', 1, 1)
130 C
131  sinpiy = sin(pi*y)
132  IF (sinpiy .EQ. 0.) CALL xermsg('SLATEC', 'GAMMA',
133  + 'X IS A NEGATIVE INTEGER', 4, 2)
134 C
135  gamma = -pi / (y*sinpiy*gamma)
136 C
137  RETURN
138  END