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
69 static const __float128 bp[] = {
70 1.0Q,
71 1.5Q,
72 };
73
74 /* log_2(1.5) */
75 static const __float128 dp_h[] = {
76 0.0,
77 5.8496250072115607565592654282227158546448E-1Q
78 };
79
80 /* Low part of log_2(1.5) */
81 static const __float128 dp_l[] = {
82 0.0,
83 1.0579781240112554492329533686862998106046E-16Q
84 };
85
86 static 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 */
97 static const __float128 LN[] =
98 {
99 -3.0779177200290054398792536829702930623200E1Q,
100 6.5135778082209159921251824580292116201640E1Q,
101 -4.6312921812152436921591152809994014413540E1Q,
102 1.2510208195629420304615674658258363295208E1Q,
103 -9.9266909031921425609179910128531667336670E-1Q
104 };
105 static 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 */
118 static const __float128 PN[] =
119 {
120 5.081801691915377692446852383385968225675E8Q,
121 9.360895299872484512023336636427675327355E6Q,
122 4.213701282274196030811629773097579432957E4Q,
123 5.201006511142748908655720086041570288182E1Q,
124 9.088368420359444263703202925095675982530E-3Q,
125 };
126 static const __float128 PD[] =
127 {
128 3.049081015149226615468111430031590411682E9Q,
129 1.069833887183886839966085436512368982758E8Q,
130 8.259257717868875207333991924545445705394E5Q,
131 1.872583833284143212651746812884298360922E3Q,
132 /* 1.0E0 */
133 };
134
135 static 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
powq(__float128 x,__float128 y)147 powq (__float128 x, __float128 y)
148 {
149 __float128 z, ax, z_h, z_l, p_h, p_l;
150 __float128 y1, t1, t2, r, s, 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 /* |y| is huge.
265 2^-16495 = 1/2 of smallest representable value.
266 If (1 - 1/131072)^y underflows, y > 1.4986e9 */
267 if (iy > 0x401d654b)
268 {
269 /* if (1 - 2^-113)^y underflows, y > 1.1873e38 */
270 if (iy > 0x407d654b)
271 {
272 if (ix <= 0x3ffeffff)
273 return (hy < 0) ? huge * huge : tiny * tiny;
274 if (ix >= 0x3fff0000)
275 return (hy > 0) ? huge * huge : tiny * tiny;
276 }
277 /* over/underflow if x is not close to one */
278 if (ix < 0x3ffeffff)
279 return (hy < 0) ? huge * huge : tiny * tiny;
280 if (ix > 0x3fff0000)
281 return (hy > 0) ? huge * huge : tiny * tiny;
282 }
283
284 ay = y > 0 ? y : -y;
285 if (ay < 0x1p-128)
286 y = y < 0 ? -0x1p-128 : 0x1p-128;
287
288 n = 0;
289 /* take care subnormal number */
290 if (ix < 0x00010000)
291 {
292 ax *= two113;
293 n -= 113;
294 o.value = ax;
295 ix = o.words32.w0;
296 }
297 n += ((ix) >> 16) - 0x3fff;
298 j = ix & 0x0000ffff;
299 /* determine interval */
300 ix = j | 0x3fff0000; /* normalize ix */
301 if (j <= 0x3988)
302 k = 0; /* |x|<sqrt(3/2) */
303 else if (j < 0xbb67)
304 k = 1; /* |x|<sqrt(3) */
305 else
306 {
307 k = 0;
308 n += 1;
309 ix -= 0x00010000;
310 }
311
312 o.value = ax;
313 o.words32.w0 = ix;
314 ax = o.value;
315
316 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
317 u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
318 v = one / (ax + bp[k]);
319 s = u * v;
320 s_h = s;
321
322 o.value = s_h;
323 o.words32.w3 = 0;
324 o.words32.w2 &= 0xf8000000;
325 s_h = o.value;
326 /* t_h=ax+bp[k] High */
327 t_h = ax + bp[k];
328 o.value = t_h;
329 o.words32.w3 = 0;
330 o.words32.w2 &= 0xf8000000;
331 t_h = o.value;
332 t_l = ax - (t_h - bp[k]);
333 s_l = v * ((u - s_h * t_h) - s_h * t_l);
334 /* compute log(ax) */
335 s2 = s * s;
336 u = LN[0] + s2 * (LN[1] + s2 * (LN[2] + s2 * (LN[3] + s2 * LN[4])));
337 v = LD[0] + s2 * (LD[1] + s2 * (LD[2] + s2 * (LD[3] + s2 * (LD[4] + s2))));
338 r = s2 * s2 * u / v;
339 r += s_l * (s_h + s);
340 s2 = s_h * s_h;
341 t_h = 3.0 + s2 + r;
342 o.value = t_h;
343 o.words32.w3 = 0;
344 o.words32.w2 &= 0xf8000000;
345 t_h = o.value;
346 t_l = r - ((t_h - 3.0) - s2);
347 /* u+v = s*(1+...) */
348 u = s_h * t_h;
349 v = s_l * t_h + t_l * s;
350 /* 2/(3log2)*(s+...) */
351 p_h = u + v;
352 o.value = p_h;
353 o.words32.w3 = 0;
354 o.words32.w2 &= 0xf8000000;
355 p_h = o.value;
356 p_l = v - (p_h - u);
357 z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */
358 z_l = cp_l * p_h + p_l * cp + dp_l[k];
359 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
360 t = (__float128) n;
361 t1 = (((z_h + z_l) + dp_h[k]) + t);
362 o.value = t1;
363 o.words32.w3 = 0;
364 o.words32.w2 &= 0xf8000000;
365 t1 = o.value;
366 t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
367
368 /* s (sign of result -ve**odd) = -1 else = 1 */
369 s = one;
370 if (((((uint32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
371 s = -one; /* (-ve)**(odd int) */
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 s * huge * huge; /* overflow */
389 else
390 {
391 if (p_l + ovt > z - p_h)
392 return s * 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 s * tiny * tiny; /* underflow */
401 else
402 {
403 if (p_l <= z - p_h)
404 return s * 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 z = scalbnq (z, n); /* subnormal output */
438 else
439 {
440 o.words32.w0 = j;
441 z = o.value;
442 }
443 return s * z;
444 }
445