GNU Octave  4.0.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
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
function r9gmit(A, X, ALGAP1, SGNGAM, ALX)
Definition: r9gmit.f:2
int exp(void) const
Definition: DET.h:64
subroutine algams(X, ALGAM, SGNGAM)
Definition: algams.f:2
octave_value log(void) const
Definition: ov.h:1190
function gami(A, X)
Definition: gami.f:2
subroutine xermsg(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
Definition: xermsg.f:2
function r9lgic(A, X, ALX)
Definition: r9lgic.f:2
subroutine xerclr
Definition: xerclr.f:2
subroutine xsgammainc(a, x, result)
Definition: xsgmainc.f:1
function alngam(X)
Definition: alngam.f:2
T abs(T x)
Definition: pr-output.cc:3062
octave_value sqrt(void) const
Definition: ov.h:1200