Mercurial > hg > CbC > CbC_gcc

comparison libquadmath/math/logq.c @ 68:561a7518be6b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update gcc-4.6
author Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
date Sun, 21 Aug 2011 07:07:55 +0900
parents
children 04ced10e8804
comparison
equal deleted inserted replaced
67:f6334be47118 68:561a7518be6b
1 /* logll.c
2 *
3 * Natural logarithm for 128-bit long double precision.
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * long double x, y, logl();
10 *
11 * y = logl( x );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Returns the base e (2.718...) logarithm of x.
18 *
19 * The argument is separated into its exponent and fractional
20 * parts. Use of a lookup table increases the speed of the routine.
21 * The program uses logarithms tabulated at intervals of 1/128 to
22 * cover the domain from approximately 0.7 to 1.4.
23 *
24 * On the interval [-1/128, +1/128] the logarithm of 1+x is approximated by
25 * log(1+x) = x - 0.5 x^2 + x^3 P(x) .
26 *
27 *
28 *
29 * ACCURACY:
30 *
31 * Relative error:
32 * arithmetic domain # trials peak rms
33 * IEEE 0.875, 1.125 100000 1.2e-34 4.1e-35
34 * IEEE 0.125, 8 100000 1.2e-34 4.1e-35
35 *
36 *
37 * WARNING:
38 *
39 * This program uses integer operations on bit fields of floating-point
40 * numbers. It does not work with data structures other than the
41 * structure assumed.
42 *
43 */
44
45 /* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>
46
47 This library is free software; you can redistribute it and/or
48 modify it under the terms of the GNU Lesser General Public
49 License as published by the Free Software Foundation; either
50 version 2.1 of the License, or (at your option) any later version.
51
52 This library is distributed in the hope that it will be useful,
53 but WITHOUT ANY WARRANTY; without even the implied warranty of
54 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
55 Lesser General Public License for more details.
56
57 You should have received a copy of the GNU Lesser General Public
58 License along with this library; if not, write to the Free Software
59 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
60
61 #include "quadmath-imp.h"
62
63 /* log(1+x) = x - .5 x^2 + x^3 l(x)
64 -.0078125 <= x <= +.0078125
65 peak relative error 1.2e-37 */
66 static const __float128
67 l3 = 3.333333333333333333333333333333336096926E-1Q,
68 l4 = -2.499999999999999999999999999486853077002E-1Q,
69 l5 = 1.999999999999999999999999998515277861905E-1Q,
70 l6 = -1.666666666666666666666798448356171665678E-1Q,
71 l7 = 1.428571428571428571428808945895490721564E-1Q,
72 l8 = -1.249999999999999987884655626377588149000E-1Q,
73 l9 = 1.111111111111111093947834982832456459186E-1Q,
74 l10 = -1.000000000000532974938900317952530453248E-1Q,
75 l11 = 9.090909090915566247008015301349979892689E-2Q,
76 l12 = -8.333333211818065121250921925397567745734E-2Q,
77 l13 = 7.692307559897661630807048686258659316091E-2Q,
78 l14 = -7.144242754190814657241902218399056829264E-2Q,
79 l15 = 6.668057591071739754844678883223432347481E-2Q;
80
81 /* Lookup table of ln(t) - (t-1)
82 t = 0.5 + (k+26)/128)
83 k = 0, ..., 91 */
84 static const __float128 logtbl[92] = {
85 -5.5345593589352099112142921677820359632418E-2Q,
86 -5.2108257402767124761784665198737642086148E-2Q,
87 -4.8991686870576856279407775480686721935120E-2Q,
88 -4.5993270766361228596215288742353061431071E-2Q,
89 -4.3110481649613269682442058976885699556950E-2Q,
90 -4.0340872319076331310838085093194799765520E-2Q,
91 -3.7682072451780927439219005993827431503510E-2Q,
92 -3.5131785416234343803903228503274262719586E-2Q,
93 -3.2687785249045246292687241862699949178831E-2Q,
94 -3.0347913785027239068190798397055267411813E-2Q,
95 -2.8110077931525797884641940838507561326298E-2Q,
96 -2.5972247078357715036426583294246819637618E-2Q,
97 -2.3932450635346084858612873953407168217307E-2Q,
98 -2.1988775689981395152022535153795155900240E-2Q,
99 -2.0139364778244501615441044267387667496733E-2Q,
100 -1.8382413762093794819267536615342902718324E-2Q,
101 -1.6716169807550022358923589720001638093023E-2Q,
102 -1.5138929457710992616226033183958974965355E-2Q,
103 -1.3649036795397472900424896523305726435029E-2Q,
104 -1.2244881690473465543308397998034325468152E-2Q,
105 -1.0924898127200937840689817557742469105693E-2Q,
106 -9.6875626072830301572839422532631079809328E-3Q,
107 -8.5313926245226231463436209313499745894157E-3Q,
108 -7.4549452072765973384933565912143044991706E-3Q,
109 -6.4568155251217050991200599386801665681310E-3Q,
110 -5.5356355563671005131126851708522185605193E-3Q,
111 -4.6900728132525199028885749289712348829878E-3Q,
112 -3.9188291218610470766469347968659624282519E-3Q,
113 -3.2206394539524058873423550293617843896540E-3Q,
114 -2.5942708080877805657374888909297113032132E-3Q,
115 -2.0385211375711716729239156839929281289086E-3Q,
116 -1.5522183228760777967376942769773768850872E-3Q,
117 -1.1342191863606077520036253234446621373191E-3Q,
118 -7.8340854719967065861624024730268350459991E-4Q,
119 -4.9869831458030115699628274852562992756174E-4Q,
120 -2.7902661731604211834685052867305795169688E-4Q,
121 -1.2335696813916860754951146082826952093496E-4Q,
122 -3.0677461025892873184042490943581654591817E-5Q,
123 #define ZERO logtbl[38]
124 0.0000000000000000000000000000000000000000E0Q,
125 -3.0359557945051052537099938863236321874198E-5Q,
126 -1.2081346403474584914595395755316412213151E-4Q,
127 -2.7044071846562177120083903771008342059094E-4Q,
128 -4.7834133324631162897179240322783590830326E-4Q,
129 -7.4363569786340080624467487620270965403695E-4Q,
130 -1.0654639687057968333207323853366578860679E-3Q,
131 -1.4429854811877171341298062134712230604279E-3Q,
132 -1.8753781835651574193938679595797367137975E-3Q,
133 -2.3618380914922506054347222273705859653658E-3Q,
134 -2.9015787624124743013946600163375853631299E-3Q,
135 -3.4938307889254087318399313316921940859043E-3Q,
136 -4.1378413103128673800485306215154712148146E-3Q,
137 -4.8328735414488877044289435125365629849599E-3Q,
138 -5.5782063183564351739381962360253116934243E-3Q,
139 -6.3731336597098858051938306767880719015261E-3Q,
140 -7.2169643436165454612058905294782949315193E-3Q,
141 -8.1090214990427641365934846191367315083867E-3Q,
142 -9.0486422112807274112838713105168375482480E-3Q,
143 -1.0035177140880864314674126398350812606841E-2Q,
144 -1.1067990155502102718064936259435676477423E-2Q,
145 -1.2146457974158024928196575103115488672416E-2Q,
146 -1.3269969823361415906628825374158424754308E-2Q,
147 -1.4437927104692837124388550722759686270765E-2Q,
148 -1.5649743073340777659901053944852735064621E-2Q,
149 -1.6904842527181702880599758489058031645317E-2Q,
150 -1.8202661505988007336096407340750378994209E-2Q,
151 -1.9542647000370545390701192438691126552961E-2Q,
152 -2.0924256670080119637427928803038530924742E-2Q,
153 -2.2346958571309108496179613803760727786257E-2Q,
154 -2.3810230892650362330447187267648486279460E-2Q,
155 -2.5313561699385640380910474255652501521033E-2Q,
156 -2.6856448685790244233704909690165496625399E-2Q,
157 -2.8438398935154170008519274953860128449036E-2Q,
158 -3.0058928687233090922411781058956589863039E-2Q,
159 -3.1717563112854831855692484086486099896614E-2Q,
160 -3.3413836095418743219397234253475252001090E-2Q,
161 -3.5147290019036555862676702093393332533702E-2Q,
162 -3.6917475563073933027920505457688955423688E-2Q,
163 -3.8723951502862058660874073462456610731178E-2Q,
164 -4.0566284516358241168330505467000838017425E-2Q,
165 -4.2444048996543693813649967076598766917965E-2Q,
166 -4.4356826869355401653098777649745233339196E-2Q,
167 -4.6304207416957323121106944474331029996141E-2Q,
168 -4.8285787106164123613318093945035804818364E-2Q,
169 -5.0301169421838218987124461766244507342648E-2Q,
170 -5.2349964705088137924875459464622098310997E-2Q,
171 -5.4431789996103111613753440311680967840214E-2Q,
172 -5.6546268881465384189752786409400404404794E-2Q,
173 -5.8693031345788023909329239565012647817664E-2Q,
174 -6.0871713627532018185577188079210189048340E-2Q,
175 -6.3081958078862169742820420185833800925568E-2Q,
176 -6.5323413029406789694910800219643791556918E-2Q,
177 -6.7595732653791419081537811574227049288168E-2Q
178 };
179
180 /* ln(2) = ln2a + ln2b with extended precision. */
181 static const __float128
182 ln2a = 6.93145751953125e-1Q,
183 ln2b = 1.4286068203094172321214581765680755001344E-6Q;
184
185 __float128
186 logq (__float128 x)
187 {
188 __float128 z, y, w;
189 ieee854_float128 u, t;
190 unsigned int m;
191 int k, e;
192
193 u.value = x;
194 m = u.words32.w0;
195
196 /* Check for IEEE special cases. */
197 k = m & 0x7fffffff;
198 /* log(0) = -infinity. */
199 if ((k | u.words32.w1 | u.words32.w2 | u.words32.w3) == 0)
200 {
201 return -0.5Q / ZERO;
202 }
203 /* log ( x < 0 ) = NaN */
204 if (m & 0x80000000)
205 {
206 return (x - x) / ZERO;
207 }
208 /* log (infinity or NaN) */
209 if (k >= 0x7fff0000)
210 {
211 return x + x;
212 }
213
214 /* Extract exponent and reduce domain to 0.703125 <= u < 1.40625 */
215 e = (int) (m >> 16) - (int) 0x3ffe;
216 m &= 0xffff;
217 u.words32.w0 = m | 0x3ffe0000;
218 m |= 0x10000;
219 /* Find lookup table index k from high order bits of the significand. */
220 if (m < 0x16800)
221 {
222 k = (m - 0xff00) >> 9;
223 /* t is the argument 0.5 + (k+26)/128
224 of the nearest item to u in the lookup table. */
225 t.words32.w0 = 0x3fff0000 + (k << 9);
226 t.words32.w1 = 0;
227 t.words32.w2 = 0;
228 t.words32.w3 = 0;
229 u.words32.w0 += 0x10000;
230 e -= 1;
231 k += 64;
232 }
233 else
234 {
235 k = (m - 0xfe00) >> 10;
236 t.words32.w0 = 0x3ffe0000 + (k << 10);
237 t.words32.w1 = 0;
238 t.words32.w2 = 0;
239 t.words32.w3 = 0;
240 }
241 /* On this interval the table is not used due to cancellation error. */
242 if ((x <= 1.0078125Q) && (x >= 0.9921875Q))
243 {
244 z = x - 1.0Q;
245 k = 64;
246 t.value = 1.0Q;
247 e = 0;
248 }
249 else
250 {
251 /* log(u) = log( t u/t ) = log(t) + log(u/t)
252 log(t) is tabulated in the lookup table.
253 Express log(u/t) = log(1+z), where z = u/t - 1 = (u-t)/t.
254 cf. Cody & Waite. */
255 z = (u.value - t.value) / t.value;
256 }
257 /* Series expansion of log(1+z). */
258 w = z * z;
259 y = ((((((((((((l15 * z
260 + l14) * z
261 + l13) * z
262 + l12) * z
263 + l11) * z
264 + l10) * z
265 + l9) * z
266 + l8) * z
267 + l7) * z
268 + l6) * z
269 + l5) * z
270 + l4) * z
271 + l3) * z * w;
272 y -= 0.5 * w;
273 y += e * ln2b; /* Base 2 exponent offset times ln(2). */
274 y += z;
275 y += logtbl[k-26]; /* log(t) - (t-1) */
276 y += (t.value - 1.0Q);
277 y += e * ln2a;
278 return y;
279 }