Blame view

fvn_fnlib/d9ln2r.f 8.42 KB
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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
        double precision function d9ln2r (x)
  c april 1978 edition.  w. fullerton c3, los alamos scientific lab.
  c
  c evaluate  dlog(1+x)  from 2-nd order with relative error accuracy so
  c that    dlog(1+x) = x - x**2/2 + x**3*d9ln2r(x)
  c
        double precision x, xbig, xmax, xmin, ln21cs(50), ln22cs(37),
       1  dcsevl, d1mach, dlog, dsqrt
        external d1mach, dcsevl, initds
  c
  c series for ln21       on the interval -6.25000e-01 to  0.
  c                                        with weighted error   1.82e-32
  c                                         log weighted error  31.74
  c                               significant figures required  31.00
  c                                    decimal places required  32.59
  c
        data ln21cs(  1) / +.1811196251 3478809875 8949530430 71 d+0     /
        data ln21cs(  2) / -.1562712319 2872462669 6251555410 78 d+0     /
        data ln21cs(  3) / +.2867630536 1557275209 5406271020 51 d-1     /
        data ln21cs(  4) / -.5558699655 9481398781 1577251267 81 d-2     /
        data ln21cs(  5) / +.1117897665 2299837657 3356662797 27 d-2     /
        data ln21cs(  6) / -.2308050898 2327947182 2992795857 05 d-3     /
        data ln21cs(  7) / +.4859885334 1100175874 6815580687 50 d-4     /
        data ln21cs(  8) / -.1039012738 8903210765 5142426333 38 d-4     /
        data ln21cs(  9) / +.2248456370 7390128494 6218049464 08 d-5     /
        data ln21cs( 10) / -.4914059273 9266484875 3278025970 91 d-6     /
        data ln21cs( 11) / +.1082825650 7077483336 6201529715 97 d-6     /
        data ln21cs( 12) / -.2402587276 3420701435 9766754167 19 d-7     /
        data ln21cs( 13) / +.5362460047 2708133762 9844432501 63 d-8     /
        data ln21cs( 14) / -.1202995136 2138772264 6716464243 77 d-8     /
        data ln21cs( 15) / +.2710788927 7591860785 6225516322 66 d-9     /
        data ln21cs( 16) / -.6132356261 8319010068 7967284306 90 d-10    /
        data ln21cs( 17) / +.1392085836 9159469857 4369085439 78 d-10    /
        data ln21cs( 18) / -.3169930033 0223494015 2830572608 83 d-11    /
        data ln21cs( 19) / +.7238375404 4307505335 2143261970 11 d-12    /
        data ln21cs( 20) / -.1657001718 4764411391 4988055062 68 d-12    /
        data ln21cs( 21) / +.3801842866 3117424257 3644226318 76 d-13    /
        data ln21cs( 22) / -.8741118929 6972700259 7244298991 37 d-14    /
        data ln21cs( 23) / +.2013561984 5055748302 1187510281 54 d-14    /
        data ln21cs( 24) / -.4646445640 9033907031 1020081544 77 d-15    /
        data ln21cs( 25) / +.1073928214 7018339453 4533385549 25 d-15    /
        data ln21cs( 26) / -.2485853461 9937794755 5340218339 60 d-16    /
        data ln21cs( 27) / +.5762019795 0800189813 8881426281 81 d-17    /
        data ln21cs( 28) / -.1337306376 9804394701 4021999580 50 d-17    /
        data ln21cs( 29) / +.3107465322 7331824966 5338071668 05 d-18    /
        data ln21cs( 30) / -.7228810408 3040539906 9019579176 27 d-19    /
        data ln21cs( 31) / +.1683378378 8037385103 3132581868 88 d-19    /
        data ln21cs( 32) / -.3923946331 2069958052 5193727399 25 d-20    /
        data ln21cs( 33) / +.9155146838 7536789746 3855286408 53 d-21    /
        data ln21cs( 34) / -.2137889532 1320159520 9820958010 02 d-21    /
        data ln21cs( 35) / +.4996450747 9047864699 8285645687 46 d-22    /
        data ln21cs( 36) / -.1168624063 6080170135 3608061474 13 d-22    /
        data ln21cs( 37) / +.2735312347 0391863775 6286867865 59 d-23    /
        data ln21cs( 38) / -.6406802508 4792111965 0503458815 99 d-24    /
        data ln21cs( 39) / +.1501629320 4334124162 9490719402 66 d-24    /
        data ln21cs( 40) / -.3521737241 0398479759 4971450026 66 d-25    /
        data ln21cs( 41) / +.8264390101 4814767012 4827333973 33 d-26    /
        data ln21cs( 42) / -.1940493027 5943401918 0366178986 66 d-26    /
        data ln21cs( 43) / +.4558788001 8841283562 4515884373 33 d-27    /
        data ln21cs( 44) / -.1071549208 7545202154 3786250239 99 d-27    /
        data ln21cs( 45) / +.2519940800 7927592978 0966741333 33 d-28    /
        data ln21cs( 46) / -.5928908840 0120969341 7504768000 00 d-29    /
        data ln21cs( 47) / +.1395586406 1057513058 2371532799 99 d-29    /
        data ln21cs( 48) / -.3286457881 3478583431 4366975999 99 d-30    /
        data ln21cs( 49) / +.7742496795 0478166247 2546986666 66 d-31    /
        data ln21cs( 50) / -.1824773566 7260887638 1252266666 66 d-31    /
  c
  c series for ln22       on the interval  0.          to  8.12500e-01
  c                                        with weighted error   6.10e-32
  c                                         log weighted error  31.21
  c                               significant figures required  30.32
  c                                    decimal places required  32.00
  c
        data ln22cs(  1) / -.2224253253 5020460829 8601522355 2 d+0      /
        data ln22cs(  2) / -.6104710010 8078623986 8010475576 4 d-1      /
        data ln22cs(  3) / +.7427235009 7503945905 1962975572 9 d-2      /
        data ln22cs(  4) / -.9335018261 6369705656 1277960639 7 d-3      /
        data ln22cs(  5) / +.1200499076 8726012833 5073128735 9 d-3      /
        data ln22cs(  6) / -.1570472295 2820041128 2335260824 3 d-4      /
        data ln22cs(  7) / +.2081874781 0512710960 5078359275 9 d-5      /
        data ln22cs(  8) / -.2789195577 6467136540 5721305137 5 d-6      /
        data ln22cs(  9) / +.3769355823 7601320584 2289513544 7 d-7      /
        data ln22cs( 10) / -.5130902896 5277112582 4058993800 3 d-8      /
        data ln22cs( 11) / +.7027141178 1506947382 0621821539 2 d-9      /
        data ln22cs( 12) / -.9674859550 1343423892 4397200513 7 d-10     /
        data ln22cs( 13) / +.1338104645 9248873065 8849644974 8 d-10     /
        data ln22cs( 14) / -.1858102603 5340639816 2845384659 1 d-11     /
        data ln22cs( 15) / +.2589294422 5279197493 0860012307 0 d-12     /
        data ln22cs( 16) / -.3619568316 1415886744 6602538217 2 d-13     /
        data ln22cs( 17) / +.5074037398 0166230880 0685891739 6 d-14     /
        data ln22cs( 18) / -.7131012977 0311273027 0093874892 7 d-15     /
        data ln22cs( 19) / +.1004490328 5545674818 5338678412 6 d-15     /
        data ln22cs( 20) / -.1417906532 1840257919 0440507528 5 d-16     /
        data ln22cs( 21) / +.2005297034 7433261178 9108639607 4 d-17     /
        data ln22cs( 22) / -.2840996662 3398033053 6539671756 7 d-18     /
        data ln22cs( 23) / +.4031469883 9690798995 9987866282 6 d-19     /
        data ln22cs( 24) / -.5729325241 8322073204 5549895679 9 d-20     /
        data ln22cs( 25) / +.8153488253 8900106758 4892873386 6 d-21     /
        data ln22cs( 26) / -.1161825588 5497217876 0602746879 9 d-21     /
        data ln22cs( 27) / +.1657516611 6625383436 5933977599 9 d-22     /
        data ln22cs( 28) / -.2367336704 7108051901 1401728000 0 d-23     /
        data ln22cs( 29) / +.3384670367 9755213860 7656959999 9 d-24     /
        data ln22cs( 30) / -.4843940829 2157182042 9639679999 9 d-25     /
        data ln22cs( 31) / +.6938759162 5142737186 7613866666 6 d-26     /
        data ln22cs( 32) / -.9948142607 0314365719 2379733333 3 d-27     /
        data ln22cs( 33) / +.1427440611 2116986106 3475200000 0 d-27     /
        data ln22cs( 34) / -.2049794721 8982349115 6650666666 6 d-28     /
        data ln22cs( 35) / +.2945648756 4013622228 8554666666 6 d-29     /
        data ln22cs( 36) / -.4235973185 1849570276 6933333333 3 d-30     /
        data ln22cs( 37) / +.6095532614 0038320401 0666666666 6 d-31     /
  c
        data ntln21, ntln22, xmin, xbig, xmax / 2*0, 3*0.d0 /
  c
        if (ntln21.ne.0) go to 10
        eps = d1mach(3)
        ntln21 = initds (ln21cs, 50, 0.1*eps)
        ntln22 = initds (ln22cs, 37, 0.1*eps)
  c
        xmin = -1.0d0 + dsqrt(d1mach(4))
        sqeps = sqrt (eps)
        xmax = 8.0/sqeps
        xmax = xmax - (dble(eps)*xmax**2 - 2.d0*dlog(xmax))
       1  / (2.d0*dble(eps)*xmax)
        xbig = 6.0/sqrt(sqeps)
        xbig = xbig - (dble(sqeps)*xbig**2 - 2.d0*dlog(xbig))
       1  / (2.d0*dble(sqeps)*xbig)
  c
   10   if (x.lt.(-.625d0) .or. x.gt.0.8125d0) go to 20
  c
        if (x.lt.0.0d0) d9ln2r = 0.375d0 + dcsevl (16.d0*x/5.d0+1.d0,
       1  ln21cs, ntln21)
        if (x.ge.0.0d0) d9ln2r = 0.375d0 + dcsevl (32.d0*x/13.d0-1.d0,
       1  ln22cs, ntln22)
        return
  c
   20   if (x.lt.xmin) call seteru (
       1  57hd9ln2r  answer lt half precision because x is too near -1,
       2  57, 1, 1)
        if (x.gt.xmax) call seteru (
       1  51hd9ln2r  no precision in answer because x is too big, 51,
       2  3, 2)
        if (x.gt.xbig) call seteru (
       1  53hd9ln2r  answer lt half precision because x is too big, 53,
       2  2, 1)
  c
        d9ln2r = (dlog(1.d0+x) - x*(1.d0 - 0.5d0*x)) / x**3
        return
  c
        end