Mercurial > octave-nkf
view libcruft/slatec-fn/dgamma.f @ 5100:5a92c3177fc6 before-gnuplot-split
[project @ 2004-12-27 17:20:38 by jwe]
| author | jwe |
|---|---|
| date | Mon, 27 Dec 2004 17:20:38 +0000 |
| parents | 30c606bec7a8 |
| children |
line wrap: on
line source
*DECK DGAMMA DOUBLE PRECISION FUNCTION DGAMMA (X) C***BEGIN PROLOGUE DGAMMA C***PURPOSE Compute the complete Gamma function. C***LIBRARY SLATEC (FNLIB) C***CATEGORY C7A C***TYPE DOUBLE PRECISION (GAMMA-S, DGAMMA-D, CGAMMA-C) C***KEYWORDS COMPLETE GAMMA FUNCTION, FNLIB, SPECIAL FUNCTIONS C***AUTHOR Fullerton, W., (LANL) C***DESCRIPTION C C DGAMMA(X) calculates the double precision complete Gamma function C for double precision argument X. C C Series for GAM on the interval 0. to 1.00000E+00 C with weighted error 5.79E-32 C log weighted error 31.24 C significant figures required 30.00 C decimal places required 32.05 C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH, D9LGMC, DCSEVL, DGAMLM, INITDS, XERMSG C***REVISION HISTORY (YYMMDD) C 770601 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890911 Removed unnecessary intrinsics. (WRB) C 890911 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 920618 Removed space from variable name. (RWC, WRB) C***END PROLOGUE DGAMMA DOUBLE PRECISION X, GAMCS(42), DXREL, PI, SINPIY, SQ2PIL, XMAX, 1 XMIN, Y, D9LGMC, DCSEVL, D1MACH LOGICAL FIRST C SAVE GAMCS, PI, SQ2PIL, NGAM, XMIN, XMAX, DXREL, FIRST DATA GAMCS( 1) / +.8571195590 9893314219 2006239994 2 D-2 / DATA GAMCS( 2) / +.4415381324 8410067571 9131577165 2 D-2 / DATA GAMCS( 3) / +.5685043681 5993633786 3266458878 9 D-1 / DATA GAMCS( 4) / -.4219835396 4185605010 1250018662 4 D-2 / DATA GAMCS( 5) / +.1326808181 2124602205 8400679635 2 D-2 / DATA GAMCS( 6) / -.1893024529 7988804325 2394702388 6 D-3 / DATA GAMCS( 7) / +.3606925327 4412452565 7808221722 5 D-4 / DATA GAMCS( 8) / -.6056761904 4608642184 8554829036 5 D-5 / DATA GAMCS( 9) / +.1055829546 3022833447 3182350909 3 D-5 / DATA GAMCS( 10) / -.1811967365 5423840482 9185589116 6 D-6 / DATA GAMCS( 11) / +.3117724964 7153222777 9025459316 9 D-7 / DATA GAMCS( 12) / -.5354219639 0196871408 7408102434 7 D-8 / DATA GAMCS( 13) / +.9193275519 8595889468 8778682594 0 D-9 / DATA GAMCS( 14) / -.1577941280 2883397617 6742327395 3 D-9 / DATA GAMCS( 15) / +.2707980622 9349545432 6654043308 9 D-10 / DATA GAMCS( 16) / -.4646818653 8257301440 8166105893 3 D-11 / DATA GAMCS( 17) / +.7973350192 0074196564 6076717535 9 D-12 / DATA GAMCS( 18) / -.1368078209 8309160257 9949917230 9 D-12 / DATA GAMCS( 19) / +.2347319486 5638006572 3347177168 8 D-13 / DATA GAMCS( 20) / -.4027432614 9490669327 6657053469 9 D-14 / DATA GAMCS( 21) / +.6910051747 3721009121 3833697525 7 D-15 / DATA GAMCS( 22) / -.1185584500 2219929070 5238712619 2 D-15 / DATA GAMCS( 23) / +.2034148542 4963739552 0102605193 2 D-16 / DATA GAMCS( 24) / -.3490054341 7174058492 7401294910 8 D-17 / DATA GAMCS( 25) / +.5987993856 4853055671 3505106602 6 D-18 / DATA GAMCS( 26) / -.1027378057 8722280744 9006977843 1 D-18 / DATA GAMCS( 27) / +.1762702816 0605298249 4275966074 8 D-19 / DATA GAMCS( 28) / -.3024320653 7353062609 5877211204 2 D-20 / DATA GAMCS( 29) / +.5188914660 2183978397 1783355050 6 D-21 / DATA GAMCS( 30) / -.8902770842 4565766924 4925160106 6 D-22 / DATA GAMCS( 31) / +.1527474068 4933426022 7459689130 6 D-22 / DATA GAMCS( 32) / -.2620731256 1873629002 5732833279 9 D-23 / DATA GAMCS( 33) / +.4496464047 8305386703 3104657066 6 D-24 / DATA GAMCS( 34) / -.7714712731 3368779117 0390152533 3 D-25 / DATA GAMCS( 35) / +.1323635453 1260440364 8657271466 6 D-25 / DATA GAMCS( 36) / -.2270999412 9429288167 0231381333 3 D-26 / DATA GAMCS( 37) / +.3896418998 0039914493 2081663999 9 D-27 / DATA GAMCS( 38) / -.6685198115 1259533277 9212799999 9 D-28 / DATA GAMCS( 39) / +.1146998663 1400243843 4761386666 6 D-28 / DATA GAMCS( 40) / -.1967938586 3451346772 9510399999 9 D-29 / DATA GAMCS( 41) / +.3376448816 5853380903 3489066666 6 D-30 / DATA GAMCS( 42) / -.5793070335 7821357846 2549333333 3 D-31 / DATA PI / 3.1415926535 8979323846 2643383279 50 D0 / DATA SQ2PIL / 0.9189385332 0467274178 0329736405 62 D0 / DATA FIRST /.TRUE./ C***FIRST EXECUTABLE STATEMENT DGAMMA IF (FIRST) THEN NGAM = INITDS (GAMCS, 42, 0.1*REAL(D1MACH(3)) ) C CALL DGAMLM (XMIN, XMAX) DXREL = SQRT(D1MACH(4)) ENDIF FIRST = .FALSE. C Y = ABS(X) IF (Y.GT.10.D0) GO TO 50 C C COMPUTE GAMMA(X) FOR -XBND .LE. X .LE. XBND. REDUCE INTERVAL AND FIND C GAMMA(1+Y) FOR 0.0 .LE. Y .LT. 1.0 FIRST OF ALL. C N = X IF (X.LT.0.D0) N = N - 1 Y = X - N N = N - 1 DGAMMA = 0.9375D0 + DCSEVL (2.D0*Y-1.D0, GAMCS, NGAM) IF (N.EQ.0) RETURN C IF (N.GT.0) GO TO 30 C C COMPUTE GAMMA(X) FOR X .LT. 1.0 C N = -N IF (X .EQ. 0.D0) CALL XERMSG ('SLATEC', 'DGAMMA', 'X IS 0', 4, 2) IF (X .LT. 0.0 .AND. X+N-2 .EQ. 0.D0) CALL XERMSG ('SLATEC', + 'DGAMMA', 'X IS A NEGATIVE INTEGER', 4, 2) IF (X .LT. (-0.5D0) .AND. ABS((X-AINT(X-0.5D0))/X) .LT. DXREL) + CALL XERMSG ('SLATEC', 'DGAMMA', + 'ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER', + 1, 1) C DO 20 I=1,N DGAMMA = DGAMMA/(X+I-1 ) 20 CONTINUE RETURN C C GAMMA(X) FOR X .GE. 2.0 AND X .LE. 10.0 C 30 DO 40 I=1,N DGAMMA = (Y+I) * DGAMMA 40 CONTINUE RETURN C C GAMMA(X) FOR ABS(X) .GT. 10.0. RECALL Y = ABS(X). C 50 IF (X .GT. XMAX) CALL XERMSG ('SLATEC', 'DGAMMA', + 'X SO BIG GAMMA OVERFLOWS', 3, 2) C DGAMMA = 0.D0 IF (X .LT. XMIN) CALL XERMSG ('SLATEC', 'DGAMMA', + 'X SO SMALL GAMMA UNDERFLOWS', 2, 1) IF (X.LT.XMIN) RETURN C DGAMMA = EXP ((Y-0.5D0)*LOG(Y) - Y + SQ2PIL + D9LGMC(Y) ) IF (X.GT.0.D0) RETURN C IF (ABS((X-AINT(X-0.5D0))/X) .LT. DXREL) CALL XERMSG ('SLATEC', + 'DGAMMA', + 'ANSWER LT HALF PRECISION, X TOO NEAR NEGATIVE INTEGER', 1, 1) C SINPIY = SIN (PI*Y) IF (SINPIY .EQ. 0.D0) CALL XERMSG ('SLATEC', 'DGAMMA', + 'X IS A NEGATIVE INTEGER', 4, 2) C DGAMMA = -PI/(Y*SINPIY*DGAMMA) C RETURN END