7 double precision a, x, result
14 DOUBLE PRECISION AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX,
15 $ bot, h, sga, sgngam, sqeps, t, d1mach,
d9gmit,
20 SAVE alneps, sqeps, bot, first
24 if (x .eq. 0.0d0)
then
26 if (a .eq. 0.0d0)
then
35 alneps = -
log(d1mach(3))
36 sqeps =
sqrt(d1mach(4))
41 IF (x .LT. 0.d0) CALL
xermsg(
'SLATEC',
'XGMAINC',
'X IS NEGATIVE'
44 IF (x.NE.0.d0) alx =
log(x)
46 IF (a.NE.0.d0) sga = sign(1.0d0, a)
47 ainta = aint(a + 0.5d0*sga)
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,
59 result =
exp(a*alx +
log(
d9gmit(a, x, algap1, sgngam, alx)))
62 30
IF (a.LT.x) go to 40
63 t =
d9lgit(a, x, dlngam(a+1.0d0))
66 result =
exp(a*alx + t)
69 40 alng =
d9lgic(a, x, alx)
74 IF (aeps.EQ.0.d0 .AND. ainta.LE.0.d0) go to 50
76 CALL
dlgams(a+1.0d0, algap1, sgngam)
77 t =
log(
abs(a)) + alng - algap1
78 IF (t.GT.alneps) go to 60
80 IF (t.GT.(-alneps)) h = 1.0d0 - sga * sgngam *
exp(t)
81 IF (
abs(h).GT.sqeps) go to 50
84 CALL
xermsg(
'SLATEC',
'XGMAINC',
'RESULT LT HALF PRECISION', 1,
94 60
IF (t.LT.bot) CALL
xerclr
95 result = -sga * sgngam *
exp(t)
subroutine dlgams(X, DLGAM, SGNGAM)
double precision function d9lgic(A, X, ALX)
double precision function d9gmit(A, X, ALGAP1, SGNGAM, ALX)
octave_value log(void) const
subroutine xgammainc(a, x, result)
subroutine xermsg(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
double precision function d9lgit(A, X, ALGAP1)
octave_value sqrt(void) const