1/*
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 *
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
8 * is preserved.
9 * ====================================================
10 */
11
12/* Expansions and modifications for 128-bit long double are
13 Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
14 and are incorporated herein by permission of the author. The author
15 reserves the right to distribute this material elsewhere under different
16 copying permissions. These modifications are distributed here under
17 the following terms:
18
19 This library is free software; you can redistribute it and/or
20 modify it under the terms of the GNU Lesser General Public
21 License as published by the Free Software Foundation; either
22 version 2.1 of the License, or (at your option) any later version.
23
24 This library is distributed in the hope that it will be useful,
25 but WITHOUT ANY WARRANTY; without even the implied warranty of
26 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
27 Lesser General Public License for more details.
28
29 You should have received a copy of the GNU Lesser General Public
30 License along with this library; if not, write to the Free Software
31 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
32
33/* powq(x,y) return x**y
34 *
35 * n
36 * Method: Let x = 2 * (1+f)
37 * 1. Compute and return log2(x) in two pieces:
38 * log2(x) = w1 + w2,
39 * where w1 has 113-53 = 60 bit trailing zeros.
40 * 2. Perform y*log2(x) = n+y' by simulating muti-precision
41 * arithmetic, where |y'|<=0.5.
42 * 3. Return x**y = 2**n*exp(y'*log2)
43 *
44 * Special cases:
45 * 1. (anything) ** 0 is 1
46 * 2. (anything) ** 1 is itself
47 * 3. (anything) ** NAN is NAN
48 * 4. NAN ** (anything except 0) is NAN
49 * 5. +-(|x| > 1) ** +INF is +INF
50 * 6. +-(|x| > 1) ** -INF is +0
51 * 7. +-(|x| < 1) ** +INF is +0
52 * 8. +-(|x| < 1) ** -INF is +INF
53 * 9. +-1 ** +-INF is NAN
54 * 10. +0 ** (+anything except 0, NAN) is +0
55 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0
56 * 12. +0 ** (-anything except 0, NAN) is +INF
57 * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
58 * 14. -0 ** (odd integer) = -( +0 ** (odd integer) )
59 * 15. +INF ** (+anything except 0,NAN) is +INF
60 * 16. +INF ** (-anything except 0,NAN) is +0
61 * 17. -INF ** (anything) = -0 ** (-anything)
62 * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
63 * 19. (-anything except 0 and inf) ** (non-integer) is NAN
64 *
65 */
66
67#include "quadmath-imp.h"
68
69static const __float128 bp[] = {
70 1.0Q,
71 1.5Q,
72};
73
74/* log_2(1.5) */
75static const __float128 dp_h[] = {
76 0.0,
77 5.8496250072115607565592654282227158546448E-1Q
78};
79
80/* Low part of log_2(1.5) */
81static const __float128 dp_l[] = {
82 0.0,
83 1.0579781240112554492329533686862998106046E-16Q
84};
85
86static const __float128 zero = 0.0Q,
87 one = 1.0Q,
88 two = 2.0Q,
89 two113 = 1.0384593717069655257060992658440192E34Q,
90 huge = 1.0e3000Q,
91 tiny = 1.0e-3000Q;
92
93/* 3/2 log x = 3 z + z^3 + z^3 (z^2 R(z^2))
94 z = (x-1)/(x+1)
95 1 <= x <= 1.25
96 Peak relative error 2.3e-37 */
97static const __float128 LN[] =
98{
99 -3.0779177200290054398792536829702930623200E1Q,
100 6.5135778082209159921251824580292116201640E1Q,
101 -4.6312921812152436921591152809994014413540E1Q,
102 1.2510208195629420304615674658258363295208E1Q,
103 -9.9266909031921425609179910128531667336670E-1Q
104};
105static const __float128 LD[] =
106{
107 -5.129862866715009066465422805058933131960E1Q,
108 1.452015077564081884387441590064272782044E2Q,
109 -1.524043275549860505277434040464085593165E2Q,
110 7.236063513651544224319663428634139768808E1Q,
111 -1.494198912340228235853027849917095580053E1Q
112 /* 1.0E0 */
113};
114
115/* exp(x) = 1 + x - x / (1 - 2 / (x - x^2 R(x^2)))
116 0 <= x <= 0.5
117 Peak relative error 5.7e-38 */
118static const __float128 PN[] =
119{
120 5.081801691915377692446852383385968225675E8Q,
121 9.360895299872484512023336636427675327355E6Q,
122 4.213701282274196030811629773097579432957E4Q,
123 5.201006511142748908655720086041570288182E1Q,
124 9.088368420359444263703202925095675982530E-3Q,
125};
126static const __float128 PD[] =
127{
128 3.049081015149226615468111430031590411682E9Q,
129 1.069833887183886839966085436512368982758E8Q,
130 8.259257717868875207333991924545445705394E5Q,
131 1.872583833284143212651746812884298360922E3Q,
132 /* 1.0E0 */
133};
134
135static const __float128
136 /* ln 2 */
137 lg2 = 6.9314718055994530941723212145817656807550E-1Q,
138 lg2_h = 6.9314718055994528622676398299518041312695E-1Q,
139 lg2_l = 2.3190468138462996154948554638754786504121E-17Q,
140 ovt = 8.0085662595372944372e-0017Q,
141 /* 2/(3*log(2)) */
142 cp = 9.6179669392597560490661645400126142495110E-1Q,
143 cp_h = 9.6179669392597555432899980587535537779331E-1Q,
144 cp_l = 5.0577616648125906047157785230014751039424E-17Q;
145
146__float128
147powq (__float128 x, __float128 y)
148{
149 __float128 z, ax, z_h, z_l, p_h, p_l;
150 __float128 y1, t1, t2, r, s, sgn, t, u, v, w;
151 __float128 s2, s_h, s_l, t_h, t_l, ay;
152 int32_t i, j, k, yisint, n;
153 uint32_t ix, iy;
154 int32_t hx, hy;
155 ieee854_float128 o, p, q;
156
157 p.value = x;
158 hx = p.words32.w0;
159 ix = hx & 0x7fffffff;
160
161 q.value = y;
162 hy = q.words32.w0;
163 iy = hy & 0x7fffffff;
164
165
166 /* y==zero: x**0 = 1 */
167 if ((iy | q.words32.w1 | q.words32.w2 | q.words32.w3) == 0)
168 return one;
169
170 /* 1.0**y = 1; -1.0**+-Inf = 1 */
171 if (x == one)
172 return one;
173 if (x == -1.0Q && iy == 0x7fff0000
174 && (q.words32.w1 | q.words32.w2 | q.words32.w3) == 0)
175 return one;
176
177 /* +-NaN return x+y */
178 if ((ix > 0x7fff0000)
179 || ((ix == 0x7fff0000)
180 && ((p.words32.w1 | p.words32.w2 | p.words32.w3) != 0))
181 || (iy > 0x7fff0000)
182 || ((iy == 0x7fff0000)
183 && ((q.words32.w1 | q.words32.w2 | q.words32.w3) != 0)))
184 return x + y;
185
186 /* determine if y is an odd int when x < 0
187 * yisint = 0 ... y is not an integer
188 * yisint = 1 ... y is an odd int
189 * yisint = 2 ... y is an even int
190 */
191 yisint = 0;
192 if (hx < 0)
193 {
194 if (iy >= 0x40700000) /* 2^113 */
195 yisint = 2; /* even integer y */
196 else if (iy >= 0x3fff0000) /* 1.0 */
197 {
198 if (floorq (y) == y)
199 {
200 z = 0.5 * y;
201 if (floorq (z) == z)
202 yisint = 2;
203 else
204 yisint = 1;
205 }
206 }
207 }
208
209 /* special value of y */
210 if ((q.words32.w1 | q.words32.w2 | q.words32.w3) == 0)
211 {
212 if (iy == 0x7fff0000) /* y is +-inf */
213 {
214 if (((ix - 0x3fff0000) | p.words32.w1 | p.words32.w2 | p.words32.w3)
215 == 0)
216 return y - y; /* +-1**inf is NaN */
217 else if (ix >= 0x3fff0000) /* (|x|>1)**+-inf = inf,0 */
218 return (hy >= 0) ? y : zero;
219 else /* (|x|<1)**-,+inf = inf,0 */
220 return (hy < 0) ? -y : zero;
221 }
222 if (iy == 0x3fff0000)
223 { /* y is +-1 */
224 if (hy < 0)
225 return one / x;
226 else
227 return x;
228 }
229 if (hy == 0x40000000)
230 return x * x; /* y is 2 */
231 if (hy == 0x3ffe0000)
232 { /* y is 0.5 */
233 if (hx >= 0) /* x >= +0 */
234 return sqrtq (x);
235 }
236 }
237
238 ax = fabsq (x);
239 /* special value of x */
240 if ((p.words32.w1 | p.words32.w2 | p.words32.w3) == 0)
241 {
242 if (ix == 0x7fff0000 || ix == 0 || ix == 0x3fff0000)
243 {
244 z = ax; /*x is +-0,+-inf,+-1 */
245 if (hy < 0)
246 z = one / z; /* z = (1/|x|) */
247 if (hx < 0)
248 {
249 if (((ix - 0x3fff0000) | yisint) == 0)
250 {
251 z = (z - z) / (z - z); /* (-1)**non-int is NaN */
252 }
253 else if (yisint == 1)
254 z = -z; /* (x<0)**odd = -(|x|**odd) */
255 }
256 return z;
257 }
258 }
259
260 /* (x<0)**(non-int) is NaN */
261 if (((((uint32_t) hx >> 31) - 1) | yisint) == 0)
262 return (x - x) / (x - x);
263
264 /* sgn (sign of result -ve**odd) = -1 else = 1 */
265 sgn = one;
266 if (((((uint32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
267 sgn = -one; /* (-ve)**(odd int) */
268
269 /* |y| is huge.
270 2^-16495 = 1/2 of smallest representable value.
271 If (1 - 1/131072)^y underflows, y > 1.4986e9 */
272 if (iy > 0x401d654b)
273 {
274 /* if (1 - 2^-113)^y underflows, y > 1.1873e38 */
275 if (iy > 0x407d654b)
276 {
277 if (ix <= 0x3ffeffff)
278 return (hy < 0) ? huge * huge : tiny * tiny;
279 if (ix >= 0x3fff0000)
280 return (hy > 0) ? huge * huge : tiny * tiny;
281 }
282 /* over/underflow if x is not close to one */
283 if (ix < 0x3ffeffff)
284 return (hy < 0) ? sgn * huge * huge : sgn * tiny * tiny;
285 if (ix > 0x3fff0000)
286 return (hy > 0) ? sgn * huge * huge : sgn * tiny * tiny;
287 }
288
289 ay = y > 0 ? y : -y;
290 if (ay < 0x1p-128)
291 y = y < 0 ? -0x1p-128 : 0x1p-128;
292
293 n = 0;
294 /* take care subnormal number */
295 if (ix < 0x00010000)
296 {
297 ax *= two113;
298 n -= 113;
299 o.value = ax;
300 ix = o.words32.w0;
301 }
302 n += ((ix) >> 16) - 0x3fff;
303 j = ix & 0x0000ffff;
304 /* determine interval */
305 ix = j | 0x3fff0000; /* normalize ix */
306 if (j <= 0x3988)
307 k = 0; /* |x|<sqrt(3/2) */
308 else if (j < 0xbb67)
309 k = 1; /* |x|<sqrt(3) */
310 else
311 {
312 k = 0;
313 n += 1;
314 ix -= 0x00010000;
315 }
316
317 o.value = ax;
318 o.words32.w0 = ix;
319 ax = o.value;
320
321 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
322 u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
323 v = one / (ax + bp[k]);
324 s = u * v;
325 s_h = s;
326
327 o.value = s_h;
328 o.words32.w3 = 0;
329 o.words32.w2 &= 0xf8000000;
330 s_h = o.value;
331 /* t_h=ax+bp[k] High */
332 t_h = ax + bp[k];
333 o.value = t_h;
334 o.words32.w3 = 0;
335 o.words32.w2 &= 0xf8000000;
336 t_h = o.value;
337 t_l = ax - (t_h - bp[k]);
338 s_l = v * ((u - s_h * t_h) - s_h * t_l);
339 /* compute log(ax) */
340 s2 = s * s;
341 u = LN[0] + s2 * (LN[1] + s2 * (LN[2] + s2 * (LN[3] + s2 * LN[4])));
342 v = LD[0] + s2 * (LD[1] + s2 * (LD[2] + s2 * (LD[3] + s2 * (LD[4] + s2))));
343 r = s2 * s2 * u / v;
344 r += s_l * (s_h + s);
345 s2 = s_h * s_h;
346 t_h = 3.0 + s2 + r;
347 o.value = t_h;
348 o.words32.w3 = 0;
349 o.words32.w2 &= 0xf8000000;
350 t_h = o.value;
351 t_l = r - ((t_h - 3.0) - s2);
352 /* u+v = s*(1+...) */
353 u = s_h * t_h;
354 v = s_l * t_h + t_l * s;
355 /* 2/(3log2)*(s+...) */
356 p_h = u + v;
357 o.value = p_h;
358 o.words32.w3 = 0;
359 o.words32.w2 &= 0xf8000000;
360 p_h = o.value;
361 p_l = v - (p_h - u);
362 z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */
363 z_l = cp_l * p_h + p_l * cp + dp_l[k];
364 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
365 t = (__float128) n;
366 t1 = (((z_h + z_l) + dp_h[k]) + t);
367 o.value = t1;
368 o.words32.w3 = 0;
369 o.words32.w2 &= 0xf8000000;
370 t1 = o.value;
371 t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
372
373 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
374 y1 = y;
375 o.value = y1;
376 o.words32.w3 = 0;
377 o.words32.w2 &= 0xf8000000;
378 y1 = o.value;
379 p_l = (y - y1) * t1 + y * t2;
380 p_h = y1 * t1;
381 z = p_l + p_h;
382 o.value = z;
383 j = o.words32.w0;
384 if (j >= 0x400d0000) /* z >= 16384 */
385 {
386 /* if z > 16384 */
387 if (((j - 0x400d0000) | o.words32.w1 | o.words32.w2 | o.words32.w3) != 0)
388 return sgn * huge * huge; /* overflow */
389 else
390 {
391 if (p_l + ovt > z - p_h)
392 return sgn * huge * huge; /* overflow */
393 }
394 }
395 else if ((j & 0x7fffffff) >= 0x400d01b9) /* z <= -16495 */
396 {
397 /* z < -16495 */
398 if (((j - 0xc00d01bc) | o.words32.w1 | o.words32.w2 | o.words32.w3)
399 != 0)
400 return sgn * tiny * tiny; /* underflow */
401 else
402 {
403 if (p_l <= z - p_h)
404 return sgn * tiny * tiny; /* underflow */
405 }
406 }
407 /* compute 2**(p_h+p_l) */
408 i = j & 0x7fffffff;
409 k = (i >> 16) - 0x3fff;
410 n = 0;
411 if (i > 0x3ffe0000)
412 { /* if |z| > 0.5, set n = [z+0.5] */
413 n = floorq (z + 0.5Q);
414 t = n;
415 p_h -= t;
416 }
417 t = p_l + p_h;
418 o.value = t;
419 o.words32.w3 = 0;
420 o.words32.w2 &= 0xf8000000;
421 t = o.value;
422 u = t * lg2_h;
423 v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
424 z = u + v;
425 w = v - (z - u);
426 /* exp(z) */
427 t = z * z;
428 u = PN[0] + t * (PN[1] + t * (PN[2] + t * (PN[3] + t * PN[4])));
429 v = PD[0] + t * (PD[1] + t * (PD[2] + t * (PD[3] + t)));
430 t1 = z - t * u / v;
431 r = (z * t1) / (t1 - two) - (w + z * w);
432 z = one - (r - z);
433 o.value = z;
434 j = o.words32.w0;
435 j += (n << 16);
436 if ((j >> 16) <= 0)
437 {
438 z = scalbnq (z, n); /* subnormal output */
439 __float128 force_underflow = z * z;
440 math_force_eval (force_underflow);
441 }
442 else
443 {
444 o.words32.w0 = j;
445 z = o.value;
446 }
447 return sgn * z;
448}
449