dshi.f
1.69 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
double precision function dshi (x)
c december 1980 edition. w. fullerton, bell labs.
c
c evaluate the hyperbolic sine integral
c shi = integral from 0 to x of sinh(t)/t dt.
c
double precision x, shics(10), xsml, absx,
1 d1mach, dcsevl, de1, dei, dsqrt
c
c series for shi on the interval 0.00000e+00 to 1.40625e-01
c with weighted error 7.11e-32
c log weighted error 31.15
c significant figures required 28.89
c decimal places required 31.65
c
data shi cs( 1) / 0.0078372685 6889009506 9520098431 7332d0/
data shi cs( 2) / 0.0039227664 9342345639 7269757442 7225d0/
data shi cs( 3) / 0.0000041346 7878876172 6674674790 8275d0/
data shi cs( 4) / 0.0000000024 7074803728 8274213514 5302d0/
data shi cs( 5) / 0.0000000000 0093792955 9076363045 7157d0/
data shi cs( 6) / 0.0000000000 0000024518 1701952086 7353d0/
data shi cs( 7) / 0.0000000000 0000000004 6741615525 7592d0/
data shi cs( 8) / 0.0000000000 0000000000 0006780307 2389d0/
data shi cs( 9) / 0.0000000000 0000000000 0000000773 1289d0/
data shi cs( 10) / 0.0000000000 0000000000 0000000000 0711d0/
c
data nshi, xsml /0, 0.0d0/
c
if (nshi.ne.0) go to 10
nshi = initds (shics, 10, 0.1*sngl(d1mach(3)))
xsml = dsqrt (d1mach(3))
c
10 absx = dabs(x)
c
dshi = x
if (absx.gt.xsml .and. absx.le.0.375d0) dshi = x*(1.d0 +
1 dcsevl (128.d0*x*x/9.d0-1.d0, shics, nshi))
if (absx.gt.0.375d0) dshi = 0.5d0*(dei(x) + de1(x))
c
return
end