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
xgmainc.f
Go to the documentation of this file.
1  subroutine xgammainc (a, x, result)
2 
3 c -- jwe, based on DGAMIT.
4 c
5 c -- Do a better job than dgami for large values of x.
6 
7  double precision a, x, result
8  intrinsic exp, log, sqrt, sign, aint
9  external dgami, dlngam, d9lgit, d9lgic, d9gmit
10 
11 C external dgamr
12 C DOUBLE PRECISION DGAMR
13 
14  DOUBLE PRECISION aeps, ainta, algap1, alneps, alng, alx,
15  $ bot, h, sga, sgngam, sqeps, t, d1mach, d9gmit,
17 
18  LOGICAL first
19 
20  SAVE alneps, sqeps, bot, first
21 
22  DATA first /.true./
23 
24  if (x .eq. 0.0d0) then
25 
26  if (a .eq. 0.0d0) then
27  result = 1.0d0
28  else
29  result = 0.0d0
30  endif
31 
32  else
33 
34  IF (first) THEN
35  alneps = -log(d1mach(3))
36  sqeps = sqrt(d1mach(4))
37  bot = log(d1mach(1))
38  ENDIF
39  first = .false.
40 C
41  IF (x .LT. 0.d0) CALL xermsg('SLATEC', 'XGMAINC', 'X IS NEGATIVE'
42  + , 2, 2)
43 C
44  IF (x.NE.0.d0) alx = log(x)
45  sga = 1.0d0
46  IF (a.NE.0.d0) sga = sign(1.0d0, a)
47  ainta = aint(a + 0.5d0*sga)
48  aeps = a - ainta
49 C
50 C IF (X.GT.0.D0) GO TO 20
51 C DGAMIT = 0.0D0
52 C IF (AINTA.GT.0.D0 .OR. AEPS.NE.0.D0) DGAMIT = DGAMR(A+1.0D0)
53 C RETURN
54 C
55  20 IF (x.GT.1.d0) go to 30
56  IF (a.GE.(-0.5d0) .OR. aeps.NE.0.d0) CALL dlgams(a+1.0d0, algap1,
57  1 sgngam)
58 C DGAMIT = D9GMIT (A, X, ALGAP1, SGNGAM, ALX)
59  result = exp(a*alx + log(d9gmit(a, x, algap1, sgngam, alx)))
60  RETURN
61 C
62  30 IF (a.LT.x) go to 40
63  t = d9lgit(a, x, dlngam(a+1.0d0))
64  IF (t.LT.bot) CALL xerclr
65 C DGAMIT = EXP (T)
66  result = exp(a*alx + t)
67  RETURN
68 C
69  40 alng = d9lgic(a, x, alx)
70 C
71 C EVALUATE DGAMIT IN TERMS OF LOG (DGAMIC (A, X))
72 C
73  h = 1.0d0
74  IF (aeps.EQ.0.d0 .AND. ainta.LE.0.d0) go to 50
75 C
76  CALL dlgams(a+1.0d0, algap1, sgngam)
77  t = log(abs(a)) + alng - algap1
78  IF (t.GT.alneps) go to 60
79 C
80  IF (t.GT.(-alneps)) h = 1.0d0 - sga * sgngam * exp(t)
81  IF (abs(h).GT.sqeps) go to 50
82 C
83  CALL xerclr
84  CALL xermsg('SLATEC', 'XGMAINC', 'RESULT LT HALF PRECISION', 1,
85  + 1)
86 C
87 C 50 T = -A*ALX + LOG(ABS(H))
88 C IF (T.LT.BOT) CALL XERCLR
89 C DGAMIT = SIGN (EXP(T), H)
90  50 result = h
91  RETURN
92 C
93 C 60 T = T - A*ALX
94  60 IF (t.LT.bot) CALL xerclr
95  result = -sga * sgngam * exp(t)
96  RETURN
97 
98  endif
99  return
100  end