GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
function csevl(X, CS, N)
Definition: csevl.f:3
subroutine gamlim(XMIN, XMAX)
Definition: gamlim.f:3
function gamma(X)
Definition: gamma.f:3
function inits(OS, NOS, ETA)
Definition: inits.f:3
octave_int< T > xmin(const octave_int< T > &x, const octave_int< T > &y)
octave_int< T > xmax(const octave_int< T > &x, const octave_int< T > &y)
static T abs(T x)
Definition: pr-output.cc:1678
real function r1mach(i)
Definition: r1mach.f:23
function r9lgmc(X)
Definition: r9lgmc.f:3
subroutine xermsg(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
Definition: xermsg.f:3