GNU Octave  4.2.1
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,
16  $ d9lgic, d9lgit, dlngam, dgami
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
subroutine dlgams(X, DLGAM, SGNGAM)
Definition: dlgams.f:2
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the base of natural logarithms The constant ex $e satisfies the equation log(e)
double precision function d9lgic(A, X, ALX)
Definition: d9lgic.f:2
double precision function d9gmit(A, X, ALGAP1, SGNGAM, ALX)
Definition: d9gmit.f:2
int exp(void) const
Definition: DET.h:66
subroutine xgammainc(a, x, result)
Definition: xgmainc.f:1
is false
Definition: cellfun.cc:398
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
Definition: Quad-opts.cc:233
subroutine xermsg(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
Definition: xermsg.f:2
double precision function d9lgit(A, X, ALGAP1)
Definition: d9lgit.f:2
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
subroutine xerclr
Definition: xerclr.f:2
octave_value sqrt(void) const
Definition: ov.h:1388