r9ln2r.f
4.38 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
function r9ln2r (x)
c april 1978 edition. w. fullerton c3, los alamos scientific lab.
c
c evaluate alog(1+x) from 2-nd order with relative error accuracy so
c that alog(1+x) = x - x**2/2 + x**3*r9ln2r(x)
c
real ln21cs(26), ln22cs(20)
external csevl, inits, r1mach
c
c series for ln21 on the interval -6.25000d-01 to 0.
c with weighted error 2.49e-17
c log weighted error 16.60
c significant figures required 15.87
c decimal places required 17.31
c
data ln21cs( 1) / .1811196251 3478810e0 /
data ln21cs( 2) / -.1562712319 2872463e0 /
data ln21cs( 3) / .0286763053 61557275e0 /
data ln21cs( 4) / -.0055586996 55948139e0 /
data ln21cs( 5) / .0011178976 65229983e0 /
data ln21cs( 6) / -.0002308050 89823279e0 /
data ln21cs( 7) / .0000485988 53341100e0 /
data ln21cs( 8) / -.0000103901 27388903e0 /
data ln21cs( 9) / .0000022484 56370739e0 /
data ln21cs(10) / -.0000004914 05927392e0 /
data ln21cs(11) / .0000001082 82565070e0 /
data ln21cs(12) / -.0000000240 25872763e0 /
data ln21cs(13) / .0000000053 62460047e0 /
data ln21cs(14) / -.0000000012 02995136e0 /
data ln21cs(15) / .0000000002 71078892e0 /
data ln21cs(16) / -.0000000000 61323562e0 /
data ln21cs(17) / .0000000000 13920858e0 /
data ln21cs(18) / -.0000000000 03169930e0 /
data ln21cs(19) / .0000000000 00723837e0 /
data ln21cs(20) / -.0000000000 00165700e0 /
data ln21cs(21) / .0000000000 00038018e0 /
data ln21cs(22) / -.0000000000 00008741e0 /
data ln21cs(23) / .0000000000 00002013e0 /
data ln21cs(24) / -.0000000000 00000464e0 /
data ln21cs(25) / .0000000000 00000107e0 /
data ln21cs(26) / -.0000000000 00000024e0 /
c
c series for ln22 on the interval 0. to 8.12500d-01
c with weighted error 1.42e-17
c log weighted error 16.85
c significant figures required 15.95
c decimal places required 17.50
c
data ln22cs( 1) / -.2224253253 5020461e0 /
data ln22cs( 2) / -.0610471001 08078624e0 /
data ln22cs( 3) / .0074272350 09750394e0 /
data ln22cs( 4) / -.0009335018 26163697e0 /
data ln22cs( 5) / .0001200499 07687260e0 /
data ln22cs( 6) / -.0000157047 22952820e0 /
data ln22cs( 7) / .0000020818 74781051e0 /
data ln22cs( 8) / -.0000002789 19557764e0 /
data ln22cs( 9) / .0000000376 93558237e0 /
data ln22cs(10) / -.0000000051 30902896e0 /
data ln22cs(11) / .0000000007 02714117e0 /
data ln22cs(12) / -.0000000000 96748595e0 /
data ln22cs(13) / .0000000000 13381046e0 /
data ln22cs(14) / -.0000000000 01858102e0 /
data ln22cs(15) / .0000000000 00258929e0 /
data ln22cs(16) / -.0000000000 00036195e0 /
data ln22cs(17) / .0000000000 00005074e0 /
data ln22cs(18) / -.0000000000 00000713e0 /
data ln22cs(19) / .0000000000 00000100e0 /
data ln22cs(20) / -.0000000000 00000014e0 /
c
data ntln21, ntln22, xmin, xbig, xmax / 2*0, 3*0.0 /
c
if (ntln21.ne.0) go to 10
eps = r1mach(3)
ntln21 = inits (ln21cs, 26, 0.1*eps)
ntln22 = inits (ln22cs, 20, 0.1*eps)
c
xmin = -1.0 + sqrt(r1mach(4))
sqeps = sqrt(eps)
xmax = 6.0/sqeps
xmax = xmax - (eps*xmax**2 - 2.0*alog(xmax)) / (2.0*eps*xmax)
xbig = 4.0/sqrt(sqeps)
xbig = xbig - (sqeps*xbig**2 - 2.0*alog(xbig)) / (2.*sqeps*xbig)
c
10 if (x.lt.(-0.625) .or. x.gt.0.8125) go to 20
c
if (x.lt.0.0) r9ln2r = 0.375 + csevl (16.*x/5.+1.0, ln21cs,
1 ntln21)
if (x.ge.0.0) r9ln2r = 0.375 + csevl (32.*x/13.-1.0, ln22cs,
1 ntln22)
return
c
20 if (x.lt.xmin) call seteru (
1 57hr9ln2r answer lt half precision because x is too near -1,
2 57, 1, 1)
if (x.gt.xmax) call seteru (
1 51hr9ln2r no precision in answer because x is too big, 51,3,2)
if (x.gt.xbig) call seteru (
1 53hr9ln2r answer lt half precision because x is too big, 53,
2 2, 1)
c
r9ln2r = (alog(1.0+x) - x*(1.0-0.5*x) ) / x**3
return
c
end