double precision function dpsi (x) c june 1977 edition. w. fullerton, c3, los alamos scientific lab. double precision x, psics(42), apsics(16), aux, dxrel, pi, xbig, 1 y, dint, dcot, dcsevl, d1mach, dlog, dsqrt external d1mach, dcot, dcsevl, dint, dlog, dsqrt, initds 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 data psi cs( 1) / -.3805708083 5217921520 4376776670 39 d-1 / data psi cs( 2) / +.4914153930 2938712748 2046996542 77 d+0 / data psi cs( 3) / -.5681574782 1244730242 8920647340 81 d-1 / data psi cs( 4) / +.8357821225 9143131362 7756507478 62 d-2 / data psi cs( 5) / -.1333232857 9943425998 0792741723 93 d-2 / data psi cs( 6) / +.2203132870 6930824892 8723979795 21 d-3 / data psi cs( 7) / -.3704023817 8456883592 8890869492 29 d-4 / data psi cs( 8) / +.6283793654 8549898933 6514187176 90 d-5 / data psi cs( 9) / -.1071263908 5061849855 2835417470 74 d-5 / data psi cs( 10) / +.1831283946 5484165805 7315898103 78 d-6 / data psi cs( 11) / -.3135350936 1808509869 0057797968 85 d-7 / data psi cs( 12) / +.5372808776 2007766260 4719191436 15 d-8 / data psi cs( 13) / -.9211681415 9784275717 8806326247 30 d-9 / data psi cs( 14) / +.1579812652 1481822782 2528840328 23 d-9 / data psi cs( 15) / -.2709864613 2380443065 4405894097 07 d-10 / data psi cs( 16) / +.4648722859 9096834872 9473195295 49 d-11 / data psi cs( 17) / -.7975272563 8303689726 5047977727 37 d-12 / data psi cs( 18) / +.1368272385 7476992249 2510538928 38 d-12 / data psi cs( 19) / -.2347515606 0658972717 3206779807 19 d-13 / data psi cs( 20) / +.4027630715 5603541107 9079250062 81 d-14 / data psi cs( 21) / -.6910251853 1179037846 5474229747 71 d-15 / data psi cs( 22) / +.1185604713 8863349552 9291395257 68 d-15 / data psi cs( 23) / -.2034168961 6261559308 1542104842 23 d-16 / data psi cs( 24) / +.3490074968 6463043850 3742329323 51 d-17 / data psi cs( 25) / -.5988014693 4976711003 0110813934 93 d-18 / data psi cs( 26) / +.1027380162 8080588258 3980057122 13 d-18 / data psi cs( 27) / -.1762704942 4561071368 3592601053 86 d-19 / data psi cs( 28) / +.3024322801 8156920457 4540354901 33 d-20 / data psi cs( 29) / -.5188916830 2092313774 2860888746 66 d-21 / data psi cs( 30) / +.8902773034 5845713905 0058874879 99 d-22 / data psi cs( 31) / -.1527474289 9426728392 8949719040 00 d-22 / data psi cs( 32) / +.2620731479 8962083136 3583180799 99 d-23 / data psi cs( 33) / -.4496464273 8220696772 5983880533 33 d-24 / data psi cs( 34) / +.7714712959 6345107028 9193642666 66 d-25 / data psi cs( 35) / -.1323635476 1887702968 1026389333 33 d-25 / data psi cs( 36) / +.2270999436 2408300091 2773119999 99 d-26 / data psi cs( 37) / -.3896419021 5374115954 4913919999 99 d-27 / data psi cs( 38) / +.6685198138 8855302310 6798933333 33 d-28 / data psi cs( 39) / -.1146998665 4920864872 5299199999 99 d-28 / data psi cs( 40) / +.1967938588 6541405920 5154133333 33 d-29 / data psi cs( 41) / -.3376448818 9750979801 9072000000 00 d-30 / data psi cs( 42) / +.5793070319 3214159246 6773333333 33 d-31 / 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 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 / c data pi / 3.1415926535 8979323846 2643383279 50 d0 / data ntpsi, ntapsi, xbig, dxrel / 2*0, 2*0.d0 / c if (ntpsi.ne.0) go to 10 ntpsi = initds (psics, 42, 0.1*sngl(d1mach(3)) ) ntapsi = initds (apsics, 16, 0.1*sngl(d1mach(3)) ) c xbig = 1.0d0/dsqrt(d1mach(3)) dxrel = dsqrt (d1mach(4)) c 10 y = dabs(x) c if (y.gt.10.0d0) go to 50 c c dpsi(x) for dabs(x) .le. 2 c n = x if (x.lt.0.d0) n = n - 1 y = x - dble(float(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 seteru (14hdpsi x is 0, 14, 2, 2) if (x.lt.0.d0 .and. x+dble(float(n-2)).eq.0.d0) call seteru ( 1 31hdpsi x is a negative integer, 31, 3, 2) if (x.lt.(-0.5d0) .and. dabs((x-dint(x-0.5d0))/x).lt.dxrel) call 1 seteru (68hdpsi answer lt half precision because x too near n 2egative integer, 68, 1, 1) c do 20 i=1,n dpsi = dpsi - 1.d0/(x+dble(float(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+dble(float(i))) 40 continue return c c dpsi(x) for dabs(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 = dlog(dabs(x)) - 0.5d0/x + aux 1 - pi*dcot(pi*x) if (x.gt.0.d0) dpsi = dlog(x) - 0.5d0/x + aux return c end