lxsd / fvn

Blame view

fvn_fnlib/cot.f 2.69 KB
edit raw normal view history
38581db0c   daniau   git-svn-id: https... 
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
 
        function cot (x)
  c may 1980 edition.   w. fullerton, c3, los alamos scientific lab.
        dimension cotcs(8)
        external  csevl, inits, r1mach
  c
  c series for cot        on the interval  0.          to  6.25000d-02
  c                                        with weighted error   3.76e-17
  c                                         log weighted error  16.42
  c                               significant figures required  15.51
  c                                    decimal places required  16.88
  c
        data cot cs( 1) /    .2402591609 8295630e0 /
        data cot cs( 2) /   -.0165330316 01500228e0 /
        data cot cs( 3) /   -.0000429983 91931724e0 /
        data cot cs( 4) /   -.0000001592 83223327e0 /
        data cot cs( 5) /   -.0000000006 19109313e0 /
        data cot cs( 6) /   -.0000000000 02430197e0 /
        data cot cs( 7) /   -.0000000000 00009560e0 /
        data cot cs( 8) /   -.0000000000 00000037e0 /
  c
  c pi2rec = 2/pi - 0.625
        data pi2rec / .01161977236 75813430 e0 /
        data nterms, xmax, xsml, xmin, sqeps / 0, 4*0.0 /
  c
        if (nterms.ne.0) go to 10
        nterms = inits (cotcs, 8, 0.1*r1mach(3))
        xmax = 1.0/r1mach(4)
        xsml = sqrt (3.0*r1mach(3))
        xmin = exp ( amax1(alog(r1mach(1)), -alog(r1mach(2))) + 0.01)
        sqeps = sqrt (r1mach(4))
  c
   10   y = abs(x)
        if (abs(x).lt.xmin) call seteru (
       1  48hcot     abs(x) is zero or so small cot overflows, 48, 2, 2)
        if (y.gt.xmax) call seteru (
       1  42hcot     no precision because abs(x) is big, 42, 3, 2)
  c
  c carefully compute y * (2/pi) = (aint(y) + rem(y)) * (.625 + pi2rec)
  c = aint(.625*y) + rem(.625*y) + y*pi2rec  =  aint(.625*y) + z
  c = aint(.625*y) + aint(z) + rem(z)
  c
        ainty = aint (y)
        yrem = y - ainty
        prodbg = 0.625*ainty
        ainty = aint (prodbg)
        y = (prodbg-ainty) + 0.625*yrem + y*pi2rec
        ainty2 = aint (y)
        ainty = ainty + ainty2
        y = y - ainty2
  c
        ifn = amod (ainty, 2.)
        if (ifn.eq.1) y = 1.0 - y
  c
        if (abs(x).gt.0.5 .and. y.lt.abs(x)*sqeps) call seteru (
       1  72hcot     answer lt half precision, abs(x) too big or x near n*
       2pi (n.ne.0), 72, 1, 1)
  c
        if (y.gt.0.25) go to 20
        if (y.eq.0.0) call seteru (29hcot     x is a multiple of pi,
       1  29, 4, 2)
        cot = 1.0/y
        if (y.gt.xsml) cot = (0.5 + csevl (32.0*y*y-1., cotcs, nterms)) /y
        go to 40
  c
   20   if (y.gt.0.5) go to 30
        cot = (0.5 + csevl (8.0*y*y-1., cotcs, nterms)) / (0.5*y)
        cot = (cot**2 - 1.0) * 0.5 / cot
        go to 40
  c
   30   cot = (0.5 + csevl (2.0*y*y-1., cotcs, nterms)) / (0.25*y)
        cot = (cot**2 - 1.0) * 0.5 / cot
        cot = (cot**2 - 1.0) * 0.5 / cot
  c
   40   if (x.ne.0.) cot = sign (cot, x)
        if (ifn.eq.1) cot = -cot
  c
        return
        end