Mercurial > hg > octave-nkf

annotate libcruft/slatec-fn/xgmainc.f @ 12747:901d466ee55a stable release-3-4-1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Version 3.4.1 released. * configure.ac (AC_INIT): Version is now 3.4.1. (OCTAVE_API_VERSION_NUMBER): Now 45. (OCTAVE_RELEASE_DATE): Now 2011-06-15.
author John W. Eaton <jwe@octave.org>
date Wed, 15 Jun 2011 10:35:37 -0400
parents c0532bbaefd7
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
4004
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
1 subroutine xgammainc (a, x, result)
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
2
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
3 c -- jwe, based on DGAMIT.
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
4 c
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
5 c -- Do a better job than dgami for large values of x.
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
6
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
7 double precision a, x, result
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
8 intrinsic exp, log, sqrt, sign, aint
4079
c0532bbaefd7 [project @ 2002-09-30 19:08:09 by jwe]
jwe
parents: 4004
diff changeset
9 external dgami, dlngam, d9lgit, d9lgic, d9gmit
c0532bbaefd7 [project @ 2002-09-30 19:08:09 by jwe]
jwe
parents: 4004
diff changeset
10
c0532bbaefd7 [project @ 2002-09-30 19:08:09 by jwe]
jwe
parents: 4004
diff changeset
11 C external dgamr
c0532bbaefd7 [project @ 2002-09-30 19:08:09 by jwe]
jwe
parents: 4004
diff changeset
12 C DOUBLE PRECISION DGAMR
4004
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
13
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
14 DOUBLE PRECISION AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX,
4079
c0532bbaefd7 [project @ 2002-09-30 19:08:09 by jwe]
jwe
parents: 4004
diff changeset
15 $ BOT, H, SGA, SGNGAM, SQEPS, T, D1MACH, D9GMIT,
4004
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
16 $ D9LGIC, D9LGIT, DLNGAM, DGAMI
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
17
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
18 LOGICAL FIRST
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
19
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
20 SAVE ALNEPS, SQEPS, BOT, FIRST
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
21
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
22 DATA FIRST /.TRUE./
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
23
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
24 if (x .eq. 0.0d0) then
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
25
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
26 if (a .eq. 0.0d0) then
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
27 result = 1.0d0
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
28 else
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
29 result = 0.0d0
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
30 endif
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
31
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
32 else
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
33
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
34 IF (FIRST) THEN
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
35 ALNEPS = -LOG (D1MACH(3))
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
36 SQEPS = SQRT(D1MACH(4))
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
37 BOT = LOG (D1MACH(1))
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
38 ENDIF
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
39 FIRST = .FALSE.
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
40 C
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
41 IF (X .LT. 0.D0) CALL XERMSG ('SLATEC', 'XGMAINC', 'X IS NEGATIVE'
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
42 + , 2, 2)
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
43 C
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
44 IF (X.NE.0.D0) ALX = LOG (X)
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
45 SGA = 1.0D0
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
46 IF (A.NE.0.D0) SGA = SIGN (1.0D0, A)
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
47 AINTA = AINT (A + 0.5D0*SGA)
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
48 AEPS = A - AINTA
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
49 C
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
50 C IF (X.GT.0.D0) GO TO 20
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
51 C DGAMIT = 0.0D0
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
52 C IF (AINTA.GT.0.D0 .OR. AEPS.NE.0.D0) DGAMIT = DGAMR(A+1.0D0)
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
53 C RETURN
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
54 C
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
55 20 IF (X.GT.1.D0) GO TO 30
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
56 IF (A.GE.(-0.5D0) .OR. AEPS.NE.0.D0) CALL DLGAMS (A+1.0D0, ALGAP1,
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
57 1 SGNGAM)
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
58 C DGAMIT = D9GMIT (A, X, ALGAP1, SGNGAM, ALX)
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
59 result = exp (a*alx + log (D9GMIT (A, X, ALGAP1, SGNGAM, ALX)))
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
60 RETURN
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
61 C
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
62 30 IF (A.LT.X) GO TO 40
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
63 T = D9LGIT (A, X, DLNGAM(A+1.0D0))
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
64 IF (T.LT.BOT) CALL XERCLR
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
65 C DGAMIT = EXP (T)
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
66 result = EXP (a*alx + T)
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
67 RETURN
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
68 C
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
69 40 ALNG = D9LGIC (A, X, ALX)
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
70 C
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
71 C EVALUATE DGAMIT IN TERMS OF LOG (DGAMIC (A, X))
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
72 C
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
73 H = 1.0D0
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
74 IF (AEPS.EQ.0.D0 .AND. AINTA.LE.0.D0) GO TO 50
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
75 C
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
76 CALL DLGAMS (A+1.0D0, ALGAP1, SGNGAM)
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
77 T = LOG (ABS(A)) + ALNG - ALGAP1
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
78 IF (T.GT.ALNEPS) GO TO 60
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
79 C
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
80 IF (T.GT.(-ALNEPS)) H = 1.0D0 - SGA * SGNGAM * EXP(T)
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
81 IF (ABS(H).GT.SQEPS) GO TO 50
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
82 C
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
83 CALL XERCLR
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
84 CALL XERMSG ('SLATEC', 'XGMAINC', 'RESULT LT HALF PRECISION', 1,
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
85 + 1)
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
86 C
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
87 C 50 T = -A*ALX + LOG(ABS(H))
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
88 C IF (T.LT.BOT) CALL XERCLR
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
89 C DGAMIT = SIGN (EXP(T), H)
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
90 50 result = H
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
91 RETURN
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
92 C
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
93 C 60 T = T - A*ALX
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
94 60 IF (T.LT.BOT) CALL XERCLR
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
95 result = -SGA * SGNGAM * EXP(T)
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
96 RETURN
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
97
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
98 endif
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
99 return
ca854fb51a88 [project @ 2002-07-25 06:31:33 by jwe]
jwe
parents:
diff changeset
100 end