2 DOUBLE PRECISION FUNCTION dgamma (X)
32 DOUBLE PRECISION X, GAMCS(42), DXREL, PI, SINPIY, SQ2PIL, XMAX,
33 1 xmin, y,
d9lgmc, dcsevl, d1mach
36 SAVE gamcs, pi, sq2pil, ngam, xmin, xmax, dxrel, first
37 DATA gamcs( 1) / +.8571195590 9893314219 2006239994 2 d-2 /
38 DATA gamcs( 2) / +.4415381324 8410067571 9131577165 2 d-2 /
39 DATA gamcs( 3) / +.5685043681 5993633786 3266458878 9 d-1 /
40 DATA gamcs( 4) / -.4219835396 4185605010 1250018662 4 d-2 /
41 DATA gamcs( 5) / +.1326808181 2124602205 8400679635 2 d-2 /
42 DATA gamcs( 6) / -.1893024529 7988804325 2394702388 6 d-3 /
43 DATA gamcs( 7) / +.3606925327 4412452565 7808221722 5 d-4 /
44 DATA gamcs( 8) / -.6056761904 4608642184 8554829036 5 d-5 /
45 DATA gamcs( 9) / +.1055829546 3022833447 3182350909 3 d-5 /
46 DATA gamcs( 10) / -.1811967365 5423840482 9185589116 6 d-6 /
47 DATA gamcs( 11) / +.3117724964 7153222777 9025459316 9 d-7 /
48 DATA gamcs( 12) / -.5354219639 0196871408 7408102434 7 d-8 /
49 DATA gamcs( 13) / +.9193275519 8595889468 8778682594 0 d-9 /
50 DATA gamcs( 14) / -.1577941280 2883397617 6742327395 3 d-9 /
51 DATA gamcs( 15) / +.2707980622 9349545432 6654043308 9 d-10 /
52 DATA gamcs( 16) / -.4646818653 8257301440 8166105893 3 d-11 /
53 DATA gamcs( 17) / +.7973350192 0074196564 6076717535 9 d-12 /
54 DATA gamcs( 18) / -.1368078209 8309160257 9949917230 9 d-12 /
55 DATA gamcs( 19) / +.2347319486 5638006572 3347177168 8 d-13 /
56 DATA gamcs( 20) / -.4027432614 9490669327 6657053469 9 d-14 /
57 DATA gamcs( 21) / +.6910051747 3721009121 3833697525 7 d-15 /
58 DATA gamcs( 22) / -.1185584500 2219929070 5238712619 2 d-15 /
59 DATA gamcs( 23) / +.2034148542 4963739552 0102605193 2 d-16 /
60 DATA gamcs( 24) / -.3490054341 7174058492 7401294910 8 d-17 /
61 DATA gamcs( 25) / +.5987993856 4853055671 3505106602 6 d-18 /
62 DATA gamcs( 26) / -.1027378057 8722280744 9006977843 1 d-18 /
63 DATA gamcs( 27) / +.1762702816 0605298249 4275966074 8 d-19 /
64 DATA gamcs( 28) / -.3024320653 7353062609 5877211204 2 d-20 /
65 DATA gamcs( 29) / +.5188914660 2183978397 1783355050 6 d-21 /
66 DATA gamcs( 30) / -.8902770842 4565766924 4925160106 6 d-22 /
67 DATA gamcs( 31) / +.1527474068 4933426022 7459689130 6 d-22 /
68 DATA gamcs( 32) / -.2620731256 1873629002 5732833279 9 d-23 /
69 DATA gamcs( 33) / +.4496464047 8305386703 3104657066 6 d-24 /
70 DATA gamcs( 34) / -.7714712731 3368779117 0390152533 3 d-25 /
71 DATA gamcs( 35) / +.1323635453 1260440364 8657271466 6 d-25 /
72 DATA gamcs( 36) / -.2270999412 9429288167 0231381333 3 d-26 /
73 DATA gamcs( 37) / +.3896418998 0039914493 2081663999 9 d-27 /
74 DATA gamcs( 38) / -.6685198115 1259533277 9212799999 9 d-28 /
75 DATA gamcs( 39) / +.1146998663 1400243843 4761386666 6 d-28 /
76 DATA gamcs( 40) / -.1967938586 3451346772 9510399999 9 d-29 /
77 DATA gamcs( 41) / +.3376448816 5853380903 3489066666 6 d-30 /
78 DATA gamcs( 42) / -.5793070335 7821357846 2549333333 3 d-31 /
79 DATA pi / 3.1415926535 8979323846 2643383279 50 d0 /
80 DATA sq2pil / 0.9189385332 0467274178 0329736405 62 d0 /
84 ngam =
initds(gamcs, 42, 0.1*
REAL(D1MACH(3)) )
87 dxrel =
sqrt(d1mach(4))
92 IF (y.GT.10.d0) go
to 50
98 IF (x.LT.0.d0) n = n - 1
101 dgamma = 0.9375d0 + dcsevl(2.d0*y-1.d0, gamcs, ngam)
109 IF (x .EQ. 0.d0) CALL
xermsg(
'SLATEC',
'DGAMMA',
'X IS 0', 4, 2)
110 IF (x .LT. 0.0 .AND. x+n-2 .EQ. 0.d0) CALL
xermsg(
'SLATEC',
111 +
'DGAMMA',
'X IS A NEGATIVE INTEGER', 4, 2)
112 IF (x .LT. (-0.5d0) .AND.
abs((x-aint(x-0.5d0))/x) .LT. dxrel)
113 + CALL
xermsg(
'SLATEC',
'DGAMMA',
114 +
'ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER',
131 50
IF (x .GT. xmax) CALL
xermsg(
'SLATEC',
'DGAMMA',
132 +
'X SO BIG GAMMA OVERFLOWS', 3, 2)
135 IF (x .LT. xmin) CALL
xermsg(
'SLATEC',
'DGAMMA',
136 +
'X SO SMALL GAMMA UNDERFLOWS', 2, 1)
137 IF (x.LT.xmin)
RETURN
140 IF (x.GT.0.d0)
RETURN
142 IF (
abs((x-aint(x-0.5d0))/x) .LT. dxrel) CALL
xermsg(
'SLATEC',
144 +
'ANSWER LT HALF PRECISION, X TOO NEAR NEGATIVE INTEGER', 1, 1)
147 IF (sinpiy .EQ. 0.d0) CALL
xermsg(
'SLATEC',
'DGAMMA',
148 +
'X IS A NEGATIVE INTEGER', 4, 2)
double precision function d9lgmc(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)
subroutine dgamlm(XMIN, XMAX)
octave_value sin(void) const
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
double precision function dgamma(X)
function initds(OS, NOS, ETA)
subroutine xermsg(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
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