2 DOUBLE PRECISION FUNCTION dgamit (A, X)
51 DOUBLE PRECISION A, X, AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX,
52 1 bot, h, sga, sgngam, sqeps, t, d1mach,
dgamr,
d9gmit,
d9lgit,
55 SAVE alneps, sqeps, bot, first
59 alneps = -
log(d1mach(3))
60 sqeps =
sqrt(d1mach(4))
65 IF (x .LT. 0.d0) CALL
xermsg(
'SLATEC',
'DGAMIT',
'X IS NEGATIVE'
68 IF (x.NE.0.d0) alx =
log(x)
70 IF (a.NE.0.d0) sga = sign(1.0d0, a)
71 ainta = aint(a + 0.5d0*sga)
74 IF (x.GT.0.d0) go
to 20
76 IF (ainta.GT.0.d0 .OR. aeps.NE.0.d0) dgamit =
dgamr(a+1.0d0)
79 20
IF (x.GT.1.d0) go
to 30
80 IF (a.GE.(-0.5d0) .OR. aeps.NE.0.d0) CALL
dlgams(a+1.0d0, algap1,
82 dgamit =
d9gmit(a, x, algap1, sgngam, alx)
85 30
IF (a.LT.x) go
to 40
86 t =
d9lgit(a, x, dlngam(a+1.0d0))
91 40 alng =
d9lgic(a, x, alx)
96 IF (aeps.EQ.0.d0 .AND. ainta.LE.0.d0) go
to 50
98 CALL
dlgams(a+1.0d0, algap1, sgngam)
99 t =
log(
abs(a)) + alng - algap1
100 IF (t.GT.alneps) go
to 60
102 IF (t.GT.(-alneps)) h = 1.0d0 - sga * sgngam *
exp(t)
103 IF (
abs(h).GT.sqeps) go
to 50
106 CALL
xermsg(
'SLATEC',
'DGAMIT',
'RESULT LT HALF PRECISION', 1,
109 50 t = -a*alx +
log(
abs(h))
111 dgamit = sign(
exp(t), h)
116 dgamit = -sga * sgngam *
exp(t)
subroutine dlgams(X, DLGAM, SGNGAM)
double precision function dgamr(X)
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)
double precision function d9gmit(A, X, ALGAP1, SGNGAM, ALX)
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
subroutine xermsg(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
double precision function d9lgit(A, X, ALGAP1)
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))<
octave_value sqrt(void) const