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
xsgmainc.f
Go to the documentation of this file.
1  subroutine xsgammainc (a, x, result)
2 
3 c -- jwe, based on GAMIT.
4 c
5 c -- Do a better job than gami for large values of x.
6 
7  real a, x, result
8  intrinsic exp, log, sqrt, sign, aint
9  external gami, alngam, r9lgit, r9lgic, r9gmit
10 
11 C external gamr
12 C real GAMR
13 
14  REAL AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX,
15  $ bot, h, sga, sgngam, sqeps, t, r1mach, r9gmit,
17 
18  LOGICAL FIRST
19 
20  SAVE alneps, sqeps, bot, first
21 
22  DATA first /.true./
23 
24  if (x .eq. 0.0e0) then
25 
26  if (a .eq. 0.0e0) then
27  result = 1.0e0
28  else
29  result = 0.0e0
30  endif
31 
32  else
33 
34  IF (first) THEN
35  alneps = -log(r1mach(3))
36  sqeps = sqrt(r1mach(4))
37  bot = log(r1mach(1))
38  ENDIF
39  first = .false.
40 C
41  IF (x .LT. 0.e0) CALL xermsg('SLATEC', 'XGMAINC', 'X IS NEGATIVE'
42  + , 2, 2)
43 C
44  IF (x.NE.0.e0) alx = log(x)
45  sga = 1.0e0
46  IF (a.NE.0.e0) sga = sign(1.0e0, a)
47  ainta = aint(a + 0.5e0*sga)
48  aeps = a - ainta
49 C
50 C IF (X.GT.0.E0) GO TO 20
51 C GAMIT = 0.0E0
52 C IF (AINTA.GT.0.E0 .OR. AEPS.NE.0.E0) GAMIT = GAMR(A+1.0E0)
53 C RETURN
54 C
55  20 IF (x.GT.1.e0) go to 30
56  IF (a.GE.(-0.5e0) .OR. aeps.NE.0.e0) CALL algams(a+1.0e0, algap1,
57  1 sgngam)
58 C GAMIT = R9GMIT (A, X, ALGAP1, SGNGAM, ALX)
59  result = exp(a*alx + log(r9gmit(a, x, algap1, sgngam, alx)))
60  RETURN
61 C
62  30 IF (a.LT.x) go to 40
63  t = r9lgit(a, x, alngam(a+1.0e0))
64  IF (t.LT.bot) CALL xerclr
65 C GAMIT = EXP (T)
66  result = exp(a*alx + t)
67  RETURN
68 C
69  40 alng = r9lgic(a, x, alx)
70 C
71 C EVALUATE GAMIT IN TERMS OF LOG (DGAMIC (A, X))
72 C
73  h = 1.0e0
74  IF (aeps.EQ.0.e0 .AND. ainta.LE.0.e0) go to 50
75 C
76  CALL algams(a+1.0e0, 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.0e0 - 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 GAMIT = 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
function r9lgit(A, X, ALGAP1)
Definition: r9lgit.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)
function r9gmit(A, X, ALGAP1, SGNGAM, ALX)
Definition: r9gmit.f:2
int exp(void) const
Definition: DET.h:66
subroutine algams(X, ALGAM, SGNGAM)
Definition: algams.f:2
function gami(A, X)
Definition: gami.f:2
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
function r9lgic(A, X, ALX)
Definition: r9lgic.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
subroutine xsgammainc(a, x, result)
Definition: xsgmainc.f:1
function alngam(X)
Definition: alngam.f:2
octave_value sqrt(void) const
Definition: ov.h:1388