Mercurial > octave-nkf

annotate libcruft/slatec-fn/d9lgit.f @ 5103:e2ed74b9bfa0 after-gnuplot-split

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2004-12-28 02:43:01 by jwe]
author jwe
date Tue, 28 Dec 2004 02:43:01 +0000
parents fe6f9bd9d0e6
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
3111
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
1 *DECK D9LGIT
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
2 DOUBLE PRECISION FUNCTION D9LGIT (A, X, ALGAP1)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
3 C***BEGIN PROLOGUE D9LGIT
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
4 C***SUBSIDIARY
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
5 C***PURPOSE Compute the logarithm of Tricomi's incomplete Gamma
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
6 C function with Perron's continued fraction for large X and
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
7 C A .GE. X.
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
8 C***LIBRARY SLATEC (FNLIB)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
9 C***CATEGORY C7E
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
10 C***TYPE DOUBLE PRECISION (R9LGIT-S, D9LGIT-D)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
11 C***KEYWORDS FNLIB, INCOMPLETE GAMMA FUNCTION, LOGARITHM,
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
12 C PERRON'S CONTINUED FRACTION, SPECIAL FUNCTIONS, TRICOMI
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
13 C***AUTHOR Fullerton, W., (LANL)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
14 C***DESCRIPTION
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
15 C
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
16 C Compute the log of Tricomi's incomplete gamma function with Perron's
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
17 C continued fraction for large X and for A .GE. X.
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
18 C
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
19 C***REFERENCES (NONE)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
20 C***ROUTINES CALLED D1MACH, XERMSG
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
21 C***REVISION HISTORY (YYMMDD)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
22 C 770701 DATE WRITTEN
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
23 C 890531 Changed all specific intrinsics to generic. (WRB)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
24 C 890531 REVISION DATE from Version 3.2
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
25 C 891214 Prologue converted to Version 4.0 format. (BAB)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
26 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
27 C 900720 Routine changed from user-callable to subsidiary. (WRB)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
28 C***END PROLOGUE D9LGIT
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
29 DOUBLE PRECISION A, X, ALGAP1, AX, A1X, EPS, FK, HSTAR, P, R, S,
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
30 1 SQEPS, T, D1MACH
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
31 LOGICAL FIRST
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
32 SAVE EPS, SQEPS, FIRST
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
33 DATA FIRST /.TRUE./
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
34 C***FIRST EXECUTABLE STATEMENT D9LGIT
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
35 IF (FIRST) THEN
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
36 EPS = 0.5D0*D1MACH(3)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
37 SQEPS = SQRT(D1MACH(4))
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
38 ENDIF
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
39 FIRST = .FALSE.
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
40 C
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
41 IF (X .LE. 0.D0 .OR. A .LT. X) CALL XERMSG ('SLATEC', 'D9LGIT',
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
42 + 'X SHOULD BE GT 0.0 AND LE A', 2, 2)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
43 C
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
44 AX = A + X
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
45 A1X = AX + 1.0D0
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
46 R = 0.D0
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
47 P = 1.D0
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
48 S = P
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
49 DO 20 K=1,200
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
50 FK = K
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
51 T = (A+FK)*X*(1.D0+R)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
52 R = T/((AX+FK)*(A1X+FK)-T)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
53 P = R*P
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
54 S = S + P
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
55 IF (ABS(P).LT.EPS*S) GO TO 30
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
56 20 CONTINUE
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
57 CALL XERMSG ('SLATEC', 'D9LGIT',
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
58 + 'NO CONVERGENCE IN 200 TERMS OF CONTINUED FRACTION', 3, 2)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
59 C
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
60 30 HSTAR = 1.0D0 - X*S/A1X
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
61 IF (HSTAR .LT. SQEPS) CALL XERMSG ('SLATEC', 'D9LGIT',
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
62 + 'RESULT LESS THAN HALF PRECISION', 1, 1)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
63 C
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
64 D9LGIT = -X - ALGAP1 - LOG(HSTAR)
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
65 RETURN
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
66 C
fe6f9bd9d0e6 [project @ 1997-11-26 07:52:06 by jwe]
jwe
parents:
diff changeset
67 END