dpsi.f
6.57 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
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, 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