1/*							logq.c
2 *
3 * Natural logarithm for __float128 precision.
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * __float128 x, y, logq();
10 *
11 * y = logq( 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 */
66static const __float128
67l3 =   3.333333333333333333333333333333336096926E-1Q,
68l4 =  -2.499999999999999999999999999486853077002E-1Q,
69l5 =   1.999999999999999999999999998515277861905E-1Q,
70l6 =  -1.666666666666666666666798448356171665678E-1Q,
71l7 =   1.428571428571428571428808945895490721564E-1Q,
72l8 =  -1.249999999999999987884655626377588149000E-1Q,
73l9 =   1.111111111111111093947834982832456459186E-1Q,
74l10 = -1.000000000000532974938900317952530453248E-1Q,
75l11 =  9.090909090915566247008015301349979892689E-2Q,
76l12 = -8.333333211818065121250921925397567745734E-2Q,
77l13 =  7.692307559897661630807048686258659316091E-2Q,
78l14 = -7.144242754190814657241902218399056829264E-2Q,
79l15 =  6.668057591071739754844678883223432347481E-2Q;
80
81/* Lookup table of ln(t) - (t-1)
82    t = 0.5 + (k+26)/128)
83    k = 0, ..., 91   */
84static 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. */
181static const __float128
182  ln2a = 6.93145751953125e-1Q,
183  ln2b = 1.4286068203094172321214581765680755001344E-6Q;
184
185__float128
186logq (__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}
280