dsi.f
2.32 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
double precision function dsi (x)
c december 1980 edition, w. fullerton, bell labs.
double precision x, sics(18), pi2, xsml, absx, f, g, cosx,
1 dcsevl, d1mach, dcos, dsin, dsqrt
external d1mach, dcsevl, initds
c
c series for si on the interval 0.00000e+00 to 1.60000e+01
c with weighted error 8.58e-32
c log weighted error 31.07
c significant figures required 30.53
c decimal places required 31.69
c
data si cs( 1) / -0.1315646598 1848419289 0427517300 0457d0/
data si cs( 2) / -0.2776578526 9736018920 4828766015 7299d0/
data si cs( 3) / 0.0354414054 8666591797 4913546471 0086d0/
data si cs( 4) / -0.0025631631 4479339776 5875278836 1530d0/
data si cs( 5) / 0.0001162365 3904970092 8126492148 2985d0/
data si cs( 6) / -0.0000035904 3272416060 4267000434 7148d0/
data si cs( 7) / 0.0000000802 3421237057 1016230865 2976d0/
data si cs( 8) / -0.0000000013 5629976925 4025064993 1846d0/
data si cs( 9) / 0.0000000000 1794407215 9973677556 7759d0/
data si cs( 10) / -0.0000000000 0019083873 4308714549 0737d0/
data si cs( 11) / 0.0000000000 0000166699 8958682433 0853d0/
data si cs( 12) / -0.0000000000 0000001217 3098836850 3042d0/
data si cs( 13) / 0.0000000000 0000000007 5418186699 3865d0/
data si cs( 14) / -0.0000000000 0000000000 0401417884 2446d0/
data si cs( 15) / 0.0000000000 0000000000 0001855369 0716d0/
data si cs( 16) / -0.0000000000 0000000000 0000007516 6966d0/
data si cs( 17) / 0.0000000000 0000000000 0000000026 9113d0/
data si cs( 18) / -0.0000000000 0000000000 0000000000 0858d0/
c
data pi2 / 1.5707963267 9489661923 1321691639 75 d0 /
data nsi, xsml /0, 0.0d0/
c
if (nsi.ne.0) go to 10
nsi = initds (sics, 18, 0.1*sngl(d1mach(3)))
xsml = dsqrt(d1mach(3))
c
10 absx = dabs(x)
if (absx.gt.4.0d0) go to 20
dsi = x
if (absx.lt.xsml) return
c
dsi = x*(0.75d0 + dcsevl ((x*x-8.0d0)*.125d0, sics, nsi))
return
c
20 call d9sifg (absx, f, g)
cosx = dcos (absx)
call erroff
dsi = pi2 - f*cosx - g*dsin(x)
if (x.lt.0.0d0) dsi = -dsi
c
return
end