Mercurial > hg > CbC > CbC_gcc

comparison libquadmath/math/powq.c @ 111:04ced10e8804

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
gcc 7
author kono
date Fri, 27 Oct 2017 22:46:09 +0900
parents 561a7518be6b
children 1830386684a0
comparison
equal deleted inserted replaced
68:561a7518be6b 111:04ced10e8804
28 28
29 You should have received a copy of the GNU Lesser General Public 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 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 */ 31 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
32 32
33 /* __ieee754_powl(x,y) return x**y 33 /* powq(x,y) return x**y
34 * 34 *
35 * n 35 * n
36 * Method: Let x = 2 * (1+f) 36 * Method: Let x = 2 * (1+f)
37 * 1. Compute and return log2(x) in two pieces: 37 * 1. Compute and return log2(x) in two pieces:
38 * log2(x) = w1 + w2, 38 * log2(x) = w1 + w2,
145 145
146 __float128 146 __float128
147 powq (__float128 x, __float128 y) 147 powq (__float128 x, __float128 y)
148 { 148 {
149 __float128 z, ax, z_h, z_l, p_h, p_l; 149 __float128 z, ax, z_h, z_l, p_h, p_l;
150 __float128 y1, t1, t2, r, s, t, u, v, w; 150 __float128 y1, t1, t2, r, s, sgn, t, u, v, w;
151 __float128 s2, s_h, s_l, t_h, t_l; 151 __float128 s2, s_h, s_l, t_h, t_l, ay;
152 int32_t i, j, k, yisint, n; 152 int32_t i, j, k, yisint, n;
153 uint32_t ix, iy; 153 uint32_t ix, iy;
154 int32_t hx, hy; 154 int32_t hx, hy;
155 ieee854_float128 o, p, q; 155 ieee854_float128 o, p, q;
156 156
259 259
260 /* (x<0)**(non-int) is NaN */ 260 /* (x<0)**(non-int) is NaN */
261 if (((((uint32_t) hx >> 31) - 1) | yisint) == 0) 261 if (((((uint32_t) hx >> 31) - 1) | yisint) == 0)
262 return (x - x) / (x - x); 262 return (x - x) / (x - x);
263 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
264 /* |y| is huge. 269 /* |y| is huge.
265 2^-16495 = 1/2 of smallest representable value. 270 2^-16495 = 1/2 of smallest representable value.
266 If (1 - 1/131072)^y underflows, y > 1.4986e9 */ 271 If (1 - 1/131072)^y underflows, y > 1.4986e9 */
267 if (iy > 0x401d654b) 272 if (iy > 0x401d654b)
268 { 273 {
274 if (ix >= 0x3fff0000) 279 if (ix >= 0x3fff0000)
275 return (hy > 0) ? huge * huge : tiny * tiny; 280 return (hy > 0) ? huge * huge : tiny * tiny;
276 } 281 }
277 /* over/underflow if x is not close to one */ 282 /* over/underflow if x is not close to one */
278 if (ix < 0x3ffeffff) 283 if (ix < 0x3ffeffff)
279 return (hy < 0) ? huge * huge : tiny * tiny; 284 return (hy < 0) ? sgn * huge * huge : sgn * tiny * tiny;
280 if (ix > 0x3fff0000) 285 if (ix > 0x3fff0000)
281 return (hy > 0) ? huge * huge : tiny * tiny; 286 return (hy > 0) ? sgn * huge * huge : sgn * tiny * tiny;
282 } 287 }
288
289 ay = y > 0 ? y : -y;
290 if (ay < 0x1p-128)
291 y = y < 0 ? -0x1p-128 : 0x1p-128;
283 292
284 n = 0; 293 n = 0;
285 /* take care subnormal number */ 294 /* take care subnormal number */
286 if (ix < 0x00010000) 295 if (ix < 0x00010000)
287 { 296 {
359 o.words32.w3 = 0; 368 o.words32.w3 = 0;
360 o.words32.w2 &= 0xf8000000; 369 o.words32.w2 &= 0xf8000000;
361 t1 = o.value; 370 t1 = o.value;
362 t2 = z_l - (((t1 - t) - dp_h[k]) - z_h); 371 t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
363 372
364 /* s (sign of result -ve**odd) = -1 else = 1 */
365 s = one;
366 if (((((uint32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
367 s = -one; /* (-ve)**(odd int) */
368
369 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ 373 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
370 y1 = y; 374 y1 = y;
371 o.value = y1; 375 o.value = y1;
372 o.words32.w3 = 0; 376 o.words32.w3 = 0;
373 o.words32.w2 &= 0xf8000000; 377 o.words32.w2 &= 0xf8000000;
379 j = o.words32.w0; 383 j = o.words32.w0;
380 if (j >= 0x400d0000) /* z >= 16384 */ 384 if (j >= 0x400d0000) /* z >= 16384 */
381 { 385 {
382 /* if z > 16384 */ 386 /* if z > 16384 */
383 if (((j - 0x400d0000) | o.words32.w1 | o.words32.w2 | o.words32.w3) != 0) 387 if (((j - 0x400d0000) | o.words32.w1 | o.words32.w2 | o.words32.w3) != 0)
384 return s * huge * huge; /* overflow */ 388 return sgn * huge * huge; /* overflow */
385 else 389 else
386 { 390 {
387 if (p_l + ovt > z - p_h) 391 if (p_l + ovt > z - p_h)
388 return s * huge * huge; /* overflow */ 392 return sgn * huge * huge; /* overflow */
389 } 393 }
390 } 394 }
391 else if ((j & 0x7fffffff) >= 0x400d01b9) /* z <= -16495 */ 395 else if ((j & 0x7fffffff) >= 0x400d01b9) /* z <= -16495 */
392 { 396 {
393 /* z < -16495 */ 397 /* z < -16495 */
394 if (((j - 0xc00d01bc) | o.words32.w1 | o.words32.w2 | o.words32.w3) 398 if (((j - 0xc00d01bc) | o.words32.w1 | o.words32.w2 | o.words32.w3)
395 != 0) 399 != 0)
396 return s * tiny * tiny; /* underflow */ 400 return sgn * tiny * tiny; /* underflow */
397 else 401 else
398 { 402 {
399 if (p_l <= z - p_h) 403 if (p_l <= z - p_h)
400 return s * tiny * tiny; /* underflow */ 404 return sgn * tiny * tiny; /* underflow */
401 } 405 }
402 } 406 }
403 /* compute 2**(p_h+p_l) */ 407 /* compute 2**(p_h+p_l) */
404 i = j & 0x7fffffff; 408 i = j & 0x7fffffff;
405 k = (i >> 16) - 0x3fff; 409 k = (i >> 16) - 0x3fff;
428 z = one - (r - z); 432 z = one - (r - z);
429 o.value = z; 433 o.value = z;
430 j = o.words32.w0; 434 j = o.words32.w0;
431 j += (n << 16); 435 j += (n << 16);
432 if ((j >> 16) <= 0) 436 if ((j >> 16) <= 0)
433 z = scalbnq (z, n); /* subnormal output */ 437 {
438 z = scalbnq (z, n); /* subnormal output */
439 __float128 force_underflow = z * z;
440 math_force_eval (force_underflow);
441 }
434 else 442 else
435 { 443 {
436 o.words32.w0 = j; 444 o.words32.w0 = j;
437 z = o.value; 445 z = o.value;
438 } 446 }
439 return s * z; 447 return sgn * z;
440 } 448 }