dcinh.f
2.48 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
double precision function dcinh (x)
c december 1980 edition. w. fullerton, bell labs.
c
c evaluate the hyperbolic cin function.
c cinh(x) = integral from 0 to x of (cosh(t)-1)/t dt.
c
double precision x, cinhcs(16), eul, xsml, xmin, absx, y,
1 dchi, dcsevl, d1mach, dlog, dsqrt
external d1mach, dchi, dcsevl, initds
c
c series for cinh on the interval 0.00000e+00 to 9.00000e+00
c with weighted error 4.84e-32
c log weighted error 31.32
c significant figures required 30.21
c decimal places required 31.92
c
data cinhcs( 1) / 0.1093291636 5207344314 0742519979 5917d0/
data cinhcs( 2) / 0.0573928847 5503796764 4532342982 5108d0/
data cinhcs( 3) / 0.0028095756 9788303534 1640420894 0774d0/
data cinhcs( 4) / 0.0000828780 8407213566 5573176506 9792d0/
data cinhcs( 5) / 0.0000016278 5961739141 8557772601 8815d0/
data cinhcs( 6) / 0.0000000227 8095192558 5661985908 3591d0/
data cinhcs( 7) / 0.0000000002 3844848424 6305925728 4002d0/
data cinhcs( 8) / 0.0000000000 0193608297 8078195747 1028d0/
data cinhcs( 9) / 0.0000000000 0001254536 9832817255 9683d0/
data cinhcs( 10) / 0.0000000000 0000006636 3744949726 2300d0/
data cinhcs( 11) / 0.0000000000 0000000029 1963926359 4744d0/
data cinhcs( 12) / 0.0000000000 0000000000 1084912395 6107d0/
data cinhcs( 13) / 0.0000000000 0000000000 0003449908 0805d0/
data cinhcs( 14) / 0.0000000000 0000000000 0000009493 6664d0/
data cinhcs( 15) / 0.0000000000 0000000000 0000000022 8291d0/
data cinhcs( 16) / 0.0000000000 0000000000 0000000000 0484d0/
c
data eul / 0.5772156649 0153286060 6512090082 40 d0 /
data ncinh, xsml, xmin /0, 2*0.0d0/
c
if (ncinh.ne.0) go to 10
ncinh = initds (cinhcs, 16, 0.1*sngl(d1mach(3)))
xsml = dsqrt (d1mach(3))
xmin = 2.0d0*dsqrt(d1mach(1))
c
10 absx = dabs(x)
if (absx.gt.3.0d0) go to 20
c
dcinh = 0.0d0
if (x.ne.0.d0 .and. absx.le.xmin) call seteru (
1 39hdcinh abs(x) so small cinh underflows, 39, 1, 0)
if (absx.le.xmin) return
c
y = -1.0d0
if (absx.gt.xsml) y = x*x/9.d0 - 1.0d0
dcinh = x*x* (0.25d0 + dcsevl (y, cinhcs, ncinh))
return
c
20 dcinh = dchi(absx) - eul - dlog(absx)
c
return
end