*DECK DPSI DOUBLE PRECISION FUNCTION DPSI (X) C***BEGIN PROLOGUE DPSI C***PURPOSE Compute the Psi (or Digamma) function. C***LIBRARY SLATEC (FNLIB) C***CATEGORY C7C C***TYPE DOUBLE PRECISION (PSI-S, DPSI-D, CPSI-C) C***KEYWORDS DIGAMMA FUNCTION, FNLIB, PSI FUNCTION, SPECIAL FUNCTIONS C***AUTHOR Fullerton, W., (LANL) C***DESCRIPTION C C DPSI calculates the double precision Psi (or Digamma) function for C double precision argument X. PSI(X) is the logarithmic derivative C of the Gamma function of X. C C Series for PSI 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.93 C decimal places required 32.05 C C C Series for APSI on the interval 0. to 1.00000E-02 C with weighted error 7.75E-33 C log weighted error 32.11 C significant figures required 28.88 C decimal places required 32.71 C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH, DCOT, DCSEVL, 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 900727 Added EXTERNAL statement. (WRB) C 920618 Removed space from variable name. (RWC, WRB) C***END PROLOGUE DPSI DOUBLE PRECISION X, PSICS(42), APSICS(16), AUX, DXREL, PI, XBIG, 1 Y, DCOT, DCSEVL, D1MACH LOGICAL FIRST EXTERNAL DCOT SAVE PSICS, APSICS, PI, NTPSI, NTAPSI, XBIG, DXREL, FIRST DATA PSICS( 1) / -.3805708083 5217921520 4376776670 39 D-1 / DATA PSICS( 2) / +.4914153930 2938712748 2046996542 77 D+0 / DATA PSICS( 3) / -.5681574782 1244730242 8920647340 81 D-1 / DATA PSICS( 4) / +.8357821225 9143131362 7756507478 62 D-2 / DATA PSICS( 5) / -.1333232857 9943425998 0792741723 93 D-2 / DATA PSICS( 6) / +.2203132870 6930824892 8723979795 21 D-3 / DATA PSICS( 7) / -.3704023817 8456883592 8890869492 29 D-4 / DATA PSICS( 8) / +.6283793654 8549898933 6514187176 90 D-5 / DATA PSICS( 9) / -.1071263908 5061849855 2835417470 74 D-5 / DATA PSICS( 10) / +.1831283946 5484165805 7315898103 78 D-6 / DATA PSICS( 11) / -.3135350936 1808509869 0057797968 85 D-7 / DATA PSICS( 12) / +.5372808776 2007766260 4719191436 15 D-8 / DATA PSICS( 13) / -.9211681415 9784275717 8806326247 30 D-9 / DATA PSICS( 14) / +.1579812652 1481822782 2528840328 23 D-9 / DATA PSICS( 15) / -.2709864613 2380443065 4405894097 07 D-10 / DATA PSICS( 16) / +.4648722859 9096834872 9473195295 49 D-11 / DATA PSICS( 17) / -.7975272563 8303689726 5047977727 37 D-12 / DATA PSICS( 18) / +.1368272385 7476992249 2510538928 38 D-12 / DATA PSICS( 19) / -.2347515606 0658972717 3206779807 19 D-13 / DATA PSICS( 20) / +.4027630715 5603541107 9079250062 81 D-14 / DATA PSICS( 21) / -.6910251853 1179037846 5474229747 71 D-15 / DATA PSICS( 22) / +.1185604713 8863349552 9291395257 68 D-15 / DATA PSICS( 23) / -.2034168961 6261559308 1542104842 23 D-16 / DATA PSICS( 24) / +.3490074968 6463043850 3742329323 51 D-17 / DATA PSICS( 25) / -.5988014693 4976711003 0110813934 93 D-18 / DATA PSICS( 26) / +.1027380162 8080588258 3980057122 13 D-18 / DATA PSICS( 27) / -.1762704942 4561071368 3592601053 86 D-19 / DATA PSICS( 28) / +.3024322801 8156920457 4540354901 33 D-20 / DATA PSICS( 29) / -.5188916830 2092313774 2860888746 66 D-21 / DATA PSICS( 30) / +.8902773034 5845713905 0058874879 99 D-22 / DATA PSICS( 31) / -.1527474289 9426728392 8949719040 00 D-22 / DATA PSICS( 32) / +.2620731479 8962083136 3583180799 99 D-23 / DATA PSICS( 33) / -.4496464273 8220696772 5983880533 33 D-24 / DATA PSICS( 34) / +.7714712959 6345107028 9193642666 66 D-25 / DATA PSICS( 35) / -.1323635476 1887702968 1026389333 33 D-25 / DATA PSICS( 36) / +.2270999436 2408300091 2773119999 99 D-26 / DATA PSICS( 37) / -.3896419021 5374115954 4913919999 99 D-27 / DATA PSICS( 38) / +.6685198138 8855302310 6798933333 33 D-28 / DATA PSICS( 39) / -.1146998665 4920864872 5299199999 99 D-28 / DATA PSICS( 40) / +.1967938588 6541405920 5154133333 33 D-29 / DATA PSICS( 41) / -.3376448818 9750979801 9072000000 00 D-30 / DATA PSICS( 42) / +.5793070319 3214159246 6773333333 33 D-31 / DATA APSICS( 1) / -.8327107910 6929076017 4456932269 D-3 / DATA APSICS( 2) / -.4162518421 9273935282 1627121990 D-3 / DATA APSICS( 3) / +.1034315609 7874129117 4463193961 D-6 / DATA APSICS( 4) / -.1214681841 3590415298 7299556365 D-9 / DATA APSICS( 5) / +.3113694319 9835615552 1240278178 D-12 / DATA APSICS( 6) / -.1364613371 9317704177 6516100945 D-14 / DATA APSICS( 7) / +.9020517513 1541656513 0837974000 D-17 / DATA APSICS( 8) / -.8315429974 2159146482 9933635466 D-19 / DATA APSICS( 9) / +.1012242570 7390725418 8479482666 D-20 / DATA APSICS( 10) / -.1562702494 3562250762 0478933333 D-22 / DATA APSICS( 11) / +.2965427168 0890389613 3226666666 D-24 / DATA APSICS( 12) / -.6746868867 6570216374 1866666666 D-26 / DATA APSICS( 13) / +.1803453116 9718990421 3333333333 D-27 / DATA APSICS( 14) / -.5569016182 4598360746 6666666666 D-29 / DATA APSICS( 15) / +.1958679226 0773625173 3333333333 D-30 / DATA APSICS( 16) / -.7751958925 2333568000 0000000000 D-32 / DATA PI / 3.1415926535 8979323846 2643383279 50 D0 / DATA FIRST /.TRUE./ C***FIRST EXECUTABLE STATEMENT DPSI IF (FIRST) THEN NTPSI = INITDS (PSICS, 42, 0.1*REAL(D1MACH(3)) ) NTAPSI = INITDS (APSICS, 16, 0.1*REAL(D1MACH(3)) ) C XBIG = 1.0D0/SQRT(D1MACH(3)) DXREL = SQRT(D1MACH(4)) ENDIF FIRST = .FALSE. C Y = ABS(X) C IF (Y.GT.10.0D0) GO TO 50 C C DPSI(X) FOR ABS(X) .LE. 2 C N = X IF (X.LT.0.D0) N = N - 1 Y = X - N N = N - 1 DPSI = DCSEVL (2.D0*Y-1.D0, PSICS, NTPSI) IF (N.EQ.0) RETURN C IF (N.GT.0) GO TO 30 C N = -N IF (X .EQ. 0.D0) CALL XERMSG ('SLATEC', 'DPSI', 'X IS 0', 2, 2) IF (X .LT. 0.D0 .AND. X+N-2 .EQ. 0.D0) CALL XERMSG ('SLATEC', + 'DPSI', 'X IS A NEGATIVE INTEGER', 3, 2) IF (X .LT. (-0.5D0) .AND. ABS((X-AINT(X-0.5D0))/X) .LT. DXREL) + CALL XERMSG ('SLATEC', 'DPSI', + 'ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER', + 1, 1) C DO 20 I=1,N DPSI = DPSI - 1.D0/(X+I-1) 20 CONTINUE RETURN C C DPSI(X) FOR X .GE. 2.0 AND X .LE. 10.0 C 30 DO 40 I=1,N DPSI = DPSI + 1.0D0/(Y+I) 40 CONTINUE RETURN C C DPSI(X) FOR ABS(X) .GT. 10.0 C 50 AUX = 0.D0 IF (Y.LT.XBIG) AUX = DCSEVL (2.D0*(10.D0/Y)**2-1.D0, APSICS, 1 NTAPSI) C IF (X.LT.0.D0) DPSI = LOG(ABS(X)) - 0.5D0/X + AUX 1 - PI*DCOT(PI*X) IF (X.GT.0.D0) DPSI = LOG(X) - 0.5D0/X + AUX RETURN C END