dcin.f
2.45 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
double precision function dcin (x)
c december 1980 edition. w. fullerton, bell labs.
double precision x, cincs(18), eul, xmin, f, g, sinx, absx,
1 dcsevl, d1mach, dcos, dlog, dsin, dsqrt
external d1mach, dcsevl, initds
c
c series for cin on the interval 0.00000e+00 to 1.60000e+01
c with weighted error 2.33e-33
c log weighted error 32.63
c significant figures required 31.92
c decimal places required 33.26
c
data cin cs( 1) / 0.3707450175 0909688741 6548012285 64992d0/
data cin cs( 2) / -0.0589357489 6364446831 9568643973 63697d0/
data cin cs( 3) / 0.0053818964 2113569124 0487453262 03340d0/
data cin cs( 4) / -0.0002986005 2841962135 3195949065 63410d0/
data cin cs( 5) / 0.0000109557 2575321620 0770310544 67306d0/
data cin cs( 6) / -0.0000002840 5454877346 6304917271 87731d0/
data cin cs( 7) / 0.0000000054 6973994875 3849124578 61806d0/
data cin cs( 8) / -0.0000000000 8124187461 3181570832 77452d0/
data cin cs( 9) / 0.0000000000 0095868593 1177066090 13181d0/
data cin cs( 10) / -0.0000000000 0000920266 0043923510 31377d0/
data cin cs( 11) / 0.0000000000 0000007325 8879990178 95024d0/
data cin cs( 12) / -0.0000000000 0000000049 1437266758 42909d0/
data cin cs( 13) / 0.0000000000 0000000000 2815777467 53902d0/
data cin cs( 14) / -0.0000000000 0000000000 0013939867 88501d0/
data cin cs( 15) / 0.0000000000 0000000000 0000060224 85646d0/
data cin cs( 16) / -0.0000000000 0000000000 0000000229 04717d0/
data cin cs( 17) / 0.0000000000 0000000000 0000000000 77273d0/
data cin cs( 18) / -0.0000000000 0000000000 0000000000 00233d0/
c
data eul / 0.5772156649 0153286060 6512090082 40 d0 /
data ncin, xmin /0, 0.d0/
c
if (ncin.ne.0) go to 10
ncin = initds (cincs, 18, 0.1*sngl(d1mach(3)))
xmin = dsqrt (d1mach(1))
c
10 dcin = 0.d0
absx = dabs(x)
if (absx.ne.0.d0 .and. absx.le.xmin) call seteru (
1 38hdcin x so small that cin underflows, 38, 1, 0)
if (absx.le.xmin) return
c
if (absx.gt.4.0d0) go to 20
dcin = dcsevl((x*x-8.d0)*.125d0, cincs, ncin)*x*x
return
c
20 call d9sifg (absx, f, g)
sinx = dsin (absx)
call erroff
dcin = -f*sinx + g*dcos(absx) + dlog(absx) + eul
c
return
end