diff options
author | Alan Modra <amodra@gmail.com> | 2013-08-17 18:24:58 +0930 |
---|---|---|
committer | Alan Modra <amodra@gmail.com> | 2013-10-04 10:32:36 +0930 |
commit | 765714cafcad7e6168518c61111f07bd955a9fee (patch) | |
tree | 27c4ab59ddc279114f3d0b7064a47317cb91c3ca /sysdeps | |
parent | 4ebd120cd983c8d2ac7a234884b3ac6805d82973 (diff) | |
download | glibc-765714cafcad7e6168518c61111f07bd955a9fee.tar.gz glibc-765714cafcad7e6168518c61111f07bd955a9fee.tar.xz glibc-765714cafcad7e6168518c61111f07bd955a9fee.zip |
PowerPC floating point little-endian [3 of 15]
http://sourceware.org/ml/libc-alpha/2013-08/msg00083.html Further replacement of ieee854 macros and unions. These files also have some optimisations for comparison against 0.0L, infinity and nan. Since the ABI specifies that the high double of an IBM long double pair is the value rounded to double, a high double of 0.0 means the low double must also be 0.0. The ABI also says that infinity and nan are encoded in the high double, with the low double unspecified. This means that tests for 0.0L, +/-Infinity and +/-NaN need only check the high double. * sysdeps/ieee754/ldbl-128ibm/e_atan2l.c (__ieee754_atan2l): Rewrite all uses of ieee854 long double macros and unions. Simplify tests for long doubles that are fully specified by the high double. * sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c (__ieee754_gammal_r): Likewise. * sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c (__ieee754_ilogbl): Likewise. Remove dead code too. * sysdeps/ieee754/ldbl-128ibm/e_jnl.c (__ieee754_jnl): Likewise. (__ieee754_ynl): Likewise. * sysdeps/ieee754/ldbl-128ibm/e_log10l.c (__ieee754_log10l): Likewise. * sysdeps/ieee754/ldbl-128ibm/e_logl.c (__ieee754_logl): Likewise. * sysdeps/ieee754/ldbl-128ibm/e_powl.c (__ieee754_powl): Likewise. Remove dead code too. * sysdeps/ieee754/ldbl-128ibm/k_tanl.c (__kernel_tanl): Likewise. * sysdeps/ieee754/ldbl-128ibm/s_expm1l.c (__expm1l): Likewise. * sysdeps/ieee754/ldbl-128ibm/s_frexpl.c (__frexpl): Likewise. * sysdeps/ieee754/ldbl-128ibm/s_isinf_nsl.c (__isinf_nsl): Likewise. Simplify. * sysdeps/ieee754/ldbl-128ibm/s_isinfl.c (___isinfl): Likewise. Simplify. * sysdeps/ieee754/ldbl-128ibm/s_log1pl.c (__log1pl): Likewise. * sysdeps/ieee754/ldbl-128ibm/s_modfl.c (__modfl): Likewise. * sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c (__nextafterl): Likewise. Comment on variable precision. * sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c (__nexttoward): Likewise. * sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c (__nexttowardf): Likewise. * sysdeps/ieee754/ldbl-128ibm/s_remquol.c (__remquol): Likewise. * sysdeps/ieee754/ldbl-128ibm/s_scalblnl.c (__scalblnl): Likewise. * sysdeps/ieee754/ldbl-128ibm/s_scalbnl.c (__scalbnl): Likewise. * sysdeps/ieee754/ldbl-128ibm/s_tanhl.c (__tanhl): Likewise. * sysdeps/powerpc/fpu/libm-test-ulps: Adjust tan_towardzero ulps.
Diffstat (limited to 'sysdeps')
22 files changed, 287 insertions, 240 deletions
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_atan2l.c b/sysdeps/ieee754/ldbl-128ibm/e_atan2l.c index 3e0535561c..b625323df3 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_atan2l.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_atan2l.c @@ -56,11 +56,15 @@ __ieee754_atan2l(long double y, long double x) { long double z; int64_t k,m,hx,hy,ix,iy; - u_int64_t lx,ly; + uint64_t lx; + double xhi, xlo, yhi; - GET_LDOUBLE_WORDS64(hx,lx,x); + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); ix = hx&0x7fffffffffffffffLL; - GET_LDOUBLE_WORDS64(hy,ly,y); + yhi = ldbl_high (y); + EXTRACT_WORDS64 (hy, yhi); iy = hy&0x7fffffffffffffffLL; if(((ix)>0x7ff0000000000000LL)|| ((iy)>0x7ff0000000000000LL)) /* x or y is NaN */ @@ -70,7 +74,7 @@ __ieee754_atan2l(long double y, long double x) m = ((hy>>63)&1)|((hx>>62)&2); /* 2*sign(x)+sign(y) */ /* when y = 0 */ - if((iy|(ly&0x7fffffffffffffffLL))==0) { + if(iy==0) { switch(m) { case 0: case 1: return y; /* atan(+-0,+anything)=+-0 */ @@ -79,7 +83,7 @@ __ieee754_atan2l(long double y, long double x) } } /* when x = 0 */ - if((ix|(lx&0x7fffffffffffffff))==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; + if(ix==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; /* when x is INF */ if(ix==0x7ff0000000000000LL) { diff --git a/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c b/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c index 90d8e3f0d2..84c13de9b8 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c @@ -122,11 +122,12 @@ long double __ieee754_gammal_r (long double x, int *signgamp) { int64_t hx; - u_int64_t lx; + double xhi; - GET_LDOUBLE_WORDS64 (hx, lx, x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); - if (((hx | lx) & 0x7fffffffffffffffLL) == 0) + if ((hx & 0x7fffffffffffffffLL) == 0) { /* Return value for x == 0 is Inf with divide by zero exception. */ *signgamp = 0; diff --git a/sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c b/sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c index 55f87ed422..aeace7c977 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c @@ -31,26 +31,24 @@ static char rcsid[] = "$NetBSD: $"; int __ieee754_ilogbl(long double x) { - int64_t hx,lx; + int64_t hx; int ix; + double xhi; - GET_LDOUBLE_WORDS64(hx,lx,x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); hx &= 0x7fffffffffffffffLL; if(hx <= 0x0010000000000000LL) { - if((hx|(lx&0x7fffffffffffffffLL))==0) + if(hx==0) return FP_ILOGB0; /* ilogbl(0) = FP_ILOGB0 */ else /* subnormal x */ - if(hx==0) { - for (ix = -1043; lx>0; lx<<=1) ix -=1; - } else { - for (ix = -1022, hx<<=11; hx>0; hx<<=1) ix -=1; - } + for (ix = -1022, hx<<=11; hx>0; hx<<=1) ix -=1; return ix; } else if (hx<0x7ff0000000000000LL) return (hx>>52)-0x3ff; else if (FP_ILOGBNAN != INT_MAX) { /* ISO C99 requires ilogbl(+-Inf) == INT_MAX. */ - if (((hx^0x7ff0000000000000LL)|lx) == 0) + if (hx==0x7ff0000000000000LL) return INT_MAX; } return FP_ILOGBNAN; diff --git a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c b/sysdeps/ieee754/ldbl-128ibm/e_jnl.c index 40012e41e1..817977da5b 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_jnl.c @@ -70,26 +70,25 @@ static const long double long double __ieee754_jnl (int n, long double x) { - u_int32_t se; + uint32_t se, lx; int32_t i, ix, sgn; long double a, b, temp, di; long double z, w; - ieee854_long_double_shape_type u; + double xhi; /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x) * Thus, J(-n,x) = J(n,-x) */ - u.value = x; - se = u.parts32.w0; + xhi = ldbl_high (x); + EXTRACT_WORDS (se, lx, xhi); ix = se & 0x7fffffff; /* if J(n,NaN) is NaN */ if (ix >= 0x7ff00000) { - if ((u.parts32.w0 & 0xfffff) | u.parts32.w1 - | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3) + if (((ix - 0x7ff00000) | lx) != 0) return x + x; } @@ -298,21 +297,20 @@ strong_alias (__ieee754_jnl, __jnl_finite) long double __ieee754_ynl (int n, long double x) { - u_int32_t se; + uint32_t se, lx; int32_t i, ix; int32_t sign; long double a, b, temp; - ieee854_long_double_shape_type u; + double xhi; - u.value = x; - se = u.parts32.w0; + xhi = ldbl_high (x); + EXTRACT_WORDS (se, lx, xhi); ix = se & 0x7fffffff; /* if Y(n,NaN) is NaN */ if (ix >= 0x7ff00000) { - if ((u.parts32.w0 & 0xfffff) | u.parts32.w1 - | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3) + if (((ix - 0x7ff00000) | lx) != 0) return x + x; } if (x <= 0.0L) @@ -377,14 +375,16 @@ __ieee754_ynl (int n, long double x) a = __ieee754_y0l (x); b = __ieee754_y1l (x); /* quit if b is -inf */ - u.value = b; - se = u.parts32.w0 & 0xfff00000; + xhi = ldbl_high (b); + GET_HIGH_WORD (se, xhi); + se &= 0xfff00000; for (i = 1; i < n && se != 0xfff00000; i++) { temp = b; b = ((long double) (i + i) / x) * b - a; - u.value = b; - se = u.parts32.w0 & 0xfff00000; + xhi = ldbl_high (b); + GET_HIGH_WORD (se, xhi); + se &= 0xfff00000; a = temp; } } diff --git a/sysdeps/ieee754/ldbl-128ibm/e_log10l.c b/sysdeps/ieee754/ldbl-128ibm/e_log10l.c index fae774cea8..1a6a4a0fa3 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_log10l.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_log10l.c @@ -182,11 +182,13 @@ __ieee754_log10l (long double x) long double z; long double y; int e; - int64_t hx, lx; + int64_t hx; + double xhi; /* Test for domain */ - GET_LDOUBLE_WORDS64 (hx, lx, x); - if (((hx & 0x7fffffffffffffffLL) | (lx & 0x7fffffffffffffffLL)) == 0) + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); + if ((hx & 0x7fffffffffffffffLL) == 0) return (-1.0L / (x - x)); if (hx < 0) return (x - x) / (x - x); diff --git a/sysdeps/ieee754/ldbl-128ibm/e_logl.c b/sysdeps/ieee754/ldbl-128ibm/e_logl.c index 15b5edfab3..b7db2b9784 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_logl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_logl.c @@ -188,18 +188,20 @@ static const long double long double __ieee754_logl(long double x) { - long double z, y, w; - ieee854_long_double_shape_type u, t; + long double z, y, w, t; unsigned int m; int k, e; + double xhi; + uint32_t hx, lx; - u.value = x; - m = u.parts32.w0; + xhi = ldbl_high (x); + EXTRACT_WORDS (hx, lx, xhi); + m = hx; /* Check for IEEE special cases. */ k = m & 0x7fffffff; /* log(0) = -infinity. */ - if ((k | u.parts32.w1 | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3) == 0) + if ((k | lx) == 0) { return -0.5L / ZERO; } @@ -219,7 +221,7 @@ __ieee754_logl(long double x) { z = x - 1.0L; k = 64; - t.value = 1.0L; + t = 1.0L; e = 0; } else @@ -236,10 +238,8 @@ __ieee754_logl(long double x) k = (m - 0xff000) >> 13; /* t is the argument 0.5 + (k+26)/128 of the nearest item to u in the lookup table. */ - t.parts32.w0 = 0x3ff00000 + (k << 13); - t.parts32.w1 = 0; - t.parts32.w2 = 0; - t.parts32.w3 = 0; + INSERT_WORDS (xhi, 0x3ff00000 + (k << 13), 0); + t = xhi; w0 += 0x100000; e -= 1; k += 64; @@ -247,17 +247,15 @@ __ieee754_logl(long double x) else { k = (m - 0xfe000) >> 14; - t.parts32.w0 = 0x3fe00000 + (k << 14); - t.parts32.w1 = 0; - t.parts32.w2 = 0; - t.parts32.w3 = 0; + INSERT_WORDS (xhi, 0x3fe00000 + (k << 14), 0); + t = xhi; } - u.value = __scalbnl (u.value, ((int) ((w0 - u.parts32.w0) * 2)) >> 21); + x = __scalbnl (x, ((int) ((w0 - hx) * 2)) >> 21); /* log(u) = log( t u/t ) = log(t) + log(u/t) log(t) is tabulated in the lookup table. Express log(u/t) = log(1+z), where z = u/t - 1 = (u-t)/t. cf. Cody & Waite. */ - z = (u.value - t.value) / t.value; + z = (x - t) / t; } /* Series expansion of log(1+z). */ w = z * z; @@ -284,7 +282,7 @@ __ieee754_logl(long double x) y += e * ln2b; /* Base 2 exponent offset times ln(2). */ y += z; y += logtbl[k-26]; /* log(t) - (t-1) */ - y += (t.value - 1.0L); + y += (t - 1.0L); y += e * ln2a; return y; } diff --git a/sysdeps/ieee754/ldbl-128ibm/e_powl.c b/sysdeps/ieee754/ldbl-128ibm/e_powl.c index 8bd35d0c88..c942f2f249 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_powl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_powl.c @@ -151,37 +151,32 @@ __ieee754_powl (long double x, long double y) long double y1, t1, t2, r, s, t, u, v, w; long double s2, s_h, s_l, t_h, t_l, ay; int32_t i, j, k, yisint, n; - u_int32_t ix, iy; - int32_t hx, hy; - ieee854_long_double_shape_type o, p, q; + uint32_t ix, iy; + int32_t hx, hy, hax; + double ohi, xhi, xlo, yhi, ylo; + uint32_t lx, ly, lj; - p.value = x; - hx = p.parts32.w0; + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS (hx, lx, xhi); ix = hx & 0x7fffffff; - q.value = y; - hy = q.parts32.w0; + ldbl_unpack (y, &yhi, &ylo); + EXTRACT_WORDS (hy, ly, yhi); iy = hy & 0x7fffffff; - /* y==zero: x**0 = 1 */ - if ((iy | q.parts32.w1 | (q.parts32.w2 & 0x7fffffff) | q.parts32.w3) == 0) + if ((iy | ly) == 0) return one; /* 1.0**y = 1; -1.0**+-Inf = 1 */ if (x == one) return one; - if (x == -1.0L && iy == 0x7ff00000 - && (q.parts32.w1 | (q.parts32.w2 & 0x7fffffff) | q.parts32.w3) == 0) + if (x == -1.0L && ((iy - 0x7ff00000) | ly) == 0) return one; /* +-NaN return x+y */ - if ((ix > 0x7ff00000) - || ((ix == 0x7ff00000) - && ((p.parts32.w1 | (p.parts32.w2 & 0x7fffffff) | p.parts32.w3) != 0)) - || (iy > 0x7ff00000) - || ((iy == 0x7ff00000) - && ((q.parts32.w1 | (q.parts32.w2 & 0x7fffffff) | q.parts32.w3) != 0))) + if ((ix >= 0x7ff00000 && ((ix - 0x7ff00000) | lx) != 0) + || (iy >= 0x7ff00000 && ((iy - 0x7ff00000) | ly) != 0)) return x + y; /* determine if y is an odd int when x < 0 @@ -192,7 +187,10 @@ __ieee754_powl (long double x, long double y) yisint = 0; if (hx < 0) { - if ((q.parts32.w2 & 0x7fffffff) >= 0x43400000) /* Low part >= 2^53 */ + uint32_t low_ye; + + GET_HIGH_WORD (low_ye, ylo); + if ((low_ye & 0x7fffffff) >= 0x43400000) /* Low part >= 2^53 */ yisint = 2; /* even integer y */ else if (iy >= 0x3ff00000) /* 1.0 */ { @@ -207,42 +205,43 @@ __ieee754_powl (long double x, long double y) } } + ax = fabsl (x); + /* special value of y */ - if ((q.parts32.w1 | (q.parts32.w2 & 0x7fffffff) | q.parts32.w3) == 0) + if (ly == 0) { - if (iy == 0x7ff00000 && q.parts32.w1 == 0) /* y is +-inf */ + if (iy == 0x7ff00000) /* y is +-inf */ { - if (((ix - 0x3ff00000) | p.parts32.w1 - | (p.parts32.w2 & 0x7fffffff) | p.parts32.w3) == 0) - return y - y; /* inf**+-1 is NaN */ - else if (ix > 0x3ff00000 || fabsl (x) > 1.0L) + if (ax > one) /* (|x|>1)**+-inf = inf,0 */ return (hy >= 0) ? y : zero; else /* (|x|<1)**-,+inf = inf,0 */ return (hy < 0) ? -y : zero; } - if (iy == 0x3ff00000) - { /* y is +-1 */ - if (hy < 0) - return one / x; - else - return x; - } - if (hy == 0x40000000) - return x * x; /* y is 2 */ - if (hy == 0x3fe00000) - { /* y is 0.5 */ - if (hx >= 0) /* x >= +0 */ - return __ieee754_sqrtl (x); + if (ylo == 0.0) + { + if (iy == 0x3ff00000) + { /* y is +-1 */ + if (hy < 0) + return one / x; + else + return x; + } + if (hy == 0x40000000) + return x * x; /* y is 2 */ + if (hy == 0x3fe00000) + { /* y is 0.5 */ + if (hx >= 0) /* x >= +0 */ + return __ieee754_sqrtl (x); + } } } - ax = fabsl (x); /* special value of x */ - if ((p.parts32.w1 | (p.parts32.w2 & 0x7fffffff) | p.parts32.w3) == 0) + if (lx == 0) { - if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) + if (ix == 0x7ff00000 || ix == 0 || (ix == 0x3ff00000 && xlo == 0.0)) { z = ax; /*x is +-0,+-inf,+-1 */ if (hy < 0) @@ -294,8 +293,8 @@ __ieee754_powl (long double x, long double y) { ax *= two113; n -= 113; - o.value = ax; - ix = o.parts32.w0; + ohi = ldbl_high (ax); + GET_HIGH_WORD (ix, ohi); } n += ((ix) >> 20) - 0x3ff; j = ix & 0x000fffff; @@ -312,26 +311,19 @@ __ieee754_powl (long double x, long double y) ix -= 0x00100000; } - o.value = ax; - o.value = __scalbnl (o.value, ((int) ((ix - o.parts32.w0) * 2)) >> 21); - ax = o.value; + ohi = ldbl_high (ax); + GET_HIGH_WORD (hax, ohi); + ax = __scalbnl (ax, ((int) ((ix - hax) * 2)) >> 21); /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */ v = one / (ax + bp[k]); s = u * v; - s_h = s; + s_h = ldbl_high (s); - o.value = s_h; - o.parts32.w3 = 0; - o.parts32.w2 = 0; - s_h = o.value; /* t_h=ax+bp[k] High */ t_h = ax + bp[k]; - o.value = t_h; - o.parts32.w3 = 0; - o.parts32.w2 = 0; - t_h = o.value; + t_h = ldbl_high (t_h); t_l = ax - (t_h - bp[k]); s_l = v * ((u - s_h * t_h) - s_h * t_l); /* compute log(ax) */ @@ -342,30 +334,21 @@ __ieee754_powl (long double x, long double y) r += s_l * (s_h + s); s2 = s_h * s_h; t_h = 3.0 + s2 + r; - o.value = t_h; - o.parts32.w3 = 0; - o.parts32.w2 = 0; - t_h = o.value; + t_h = ldbl_high (t_h); t_l = r - ((t_h - 3.0) - s2); /* u+v = s*(1+...) */ u = s_h * t_h; v = s_l * t_h + t_l * s; /* 2/(3log2)*(s+...) */ p_h = u + v; - o.value = p_h; - o.parts32.w3 = 0; - o.parts32.w2 = 0; - p_h = o.value; + p_h = ldbl_high (p_h); p_l = v - (p_h - u); z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */ z_l = cp_l * p_h + p_l * cp + dp_l[k]; /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */ t = (long double) n; t1 = (((z_h + z_l) + dp_h[k]) + t); - o.value = t1; - o.parts32.w3 = 0; - o.parts32.w2 = 0; - t1 = o.value; + t1 = ldbl_high (t1); t2 = z_l - (((t1 - t) - dp_h[k]) - z_h); /* s (sign of result -ve**odd) = -1 else = 1 */ @@ -374,21 +357,16 @@ __ieee754_powl (long double x, long double y) s = -one; /* (-ve)**(odd int) */ /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ - y1 = y; - o.value = y1; - o.parts32.w3 = 0; - o.parts32.w2 = 0; - y1 = o.value; + y1 = ldbl_high (y); p_l = (y - y1) * t1 + y * t2; p_h = y1 * t1; z = p_l + p_h; - o.value = z; - j = o.parts32.w0; + ohi = ldbl_high (z); + EXTRACT_WORDS (j, lj, ohi); if (j >= 0x40d00000) /* z >= 16384 */ { /* if z > 16384 */ - if (((j - 0x40d00000) | o.parts32.w1 - | (o.parts32.w2 & 0x7fffffff) | o.parts32.w3) != 0) + if (((j - 0x40d00000) | lj) != 0) return s * huge * huge; /* overflow */ else { @@ -399,8 +377,7 @@ __ieee754_powl (long double x, long double y) else if ((j & 0x7fffffff) >= 0x40d01b90) /* z <= -16495 */ { /* z < -16495 */ - if (((j - 0xc0d01bc0) | o.parts32.w1 - | (o.parts32.w2 & 0x7fffffff) | o.parts32.w3) != 0) + if (((j - 0xc0d01bc0) | lj) != 0) return s * tiny * tiny; /* underflow */ else { @@ -419,10 +396,7 @@ __ieee754_powl (long double x, long double y) p_h -= t; } t = p_l + p_h; - o.value = t; - o.parts32.w3 = 0; - o.parts32.w2 = 0; - t = o.value; + t = ldbl_high (t); u = t * lg2_h; v = (p_l - (t - p_h)) * lg2 + t * lg2_l; z = u + v; diff --git a/sysdeps/ieee754/ldbl-128ibm/k_tanl.c b/sysdeps/ieee754/ldbl-128ibm/k_tanl.c index 1f6bad241b..bcf8b5e7d6 100644 --- a/sysdeps/ieee754/ldbl-128ibm/k_tanl.c +++ b/sysdeps/ieee754/ldbl-128ibm/k_tanl.c @@ -85,17 +85,17 @@ long double __kernel_tanl (long double x, long double y, int iy) { long double z, r, v, w, s; - int32_t ix, sign; - ieee854_long_double_shape_type u, u1; + int32_t ix, sign, hx, lx; + double xhi; - u.value = x; - ix = u.parts32.w0 & 0x7fffffff; + xhi = ldbl_high (x); + EXTRACT_WORDS (hx, lx, xhi); + ix = hx & 0x7fffffff; if (ix < 0x3c600000) /* x < 2**-57 */ { - if ((int) x == 0) - { /* generate inexact */ - if ((ix | u.parts32.w1 | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3 - | (iy + 1)) == 0) + if ((int) x == 0) /* generate inexact */ + { + if ((ix | lx | (iy + 1)) == 0) return one / fabs (x); else return (iy == 1) ? x : -one / x; @@ -103,7 +103,7 @@ __kernel_tanl (long double x, long double y, int iy) } if (ix >= 0x3fe59420) /* |x| >= 0.6743316650390625 */ { - if ((u.parts32.w0 & 0x80000000) != 0) + if ((hx & 0x80000000) != 0) { x = -x; y = -y; @@ -139,15 +139,13 @@ __kernel_tanl (long double x, long double y, int iy) { /* if allow error up to 2 ulp, simply return -1.0/(x+r) here */ /* compute -1.0/(x+r) accurately */ - u1.value = w; - u1.parts32.w2 = 0; - u1.parts32.w3 = 0; - v = r - (u1.value - x); /* u1+v = r+x */ + long double u1, z1; + + u1 = ldbl_high (w); + v = r - (u1 - x); /* u1+v = r+x */ z = -1.0 / w; - u.value = z; - u.parts32.w2 = 0; - u.parts32.w3 = 0; - s = 1.0 + u.value * u1.value; - return u.value + z * (s + u.value * v); + z1 = ldbl_high (z); + s = 1.0 + z1 * u1; + return z1 + z * (s + z1 * v); } } diff --git a/sysdeps/ieee754/ldbl-128ibm/s_expm1l.c b/sysdeps/ieee754/ldbl-128ibm/s_expm1l.c index 8808dcd896..007e785346 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_expm1l.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_expm1l.c @@ -92,19 +92,19 @@ long double __expm1l (long double x) { long double px, qx, xx; - int32_t ix, sign; - ieee854_long_double_shape_type u; + int32_t ix, lx, sign; int k; + double xhi; /* Detect infinity and NaN. */ - u.value = x; - ix = u.parts32.w0; + xhi = ldbl_high (x); + EXTRACT_WORDS (ix, lx, xhi); sign = ix & 0x80000000; ix &= 0x7fffffff; if (ix >= 0x7ff00000) { /* Infinity. */ - if (((ix & 0xfffff) | u.parts32.w1 | (u.parts32.w2&0x7fffffff) | u.parts32.w3) == 0) + if (((ix - 0x7ff00000) | lx) == 0) { if (sign) return -1.0L; @@ -116,7 +116,7 @@ __expm1l (long double x) } /* expm1(+- 0) = +- 0. */ - if ((ix == 0) && (u.parts32.w1 | (u.parts32.w2&0x7fffffff) | u.parts32.w3) == 0) + if ((ix | lx) == 0) return x; /* Overflow. */ diff --git a/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c b/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c index 3ac5374116..7e40663fd4 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c @@ -36,16 +36,21 @@ two107 = 162259276829213363391578010288128.0; /* 0x4670000000000000, 0 */ long double __frexpl(long double x, int *eptr) { - u_int64_t hx, lx, ix, ixl; + uint64_t hx, lx, ix, ixl; int64_t explo; - GET_LDOUBLE_WORDS64(hx,lx,x); + double xhi, xlo; + + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); ixl = 0x7fffffffffffffffULL&lx; ix = 0x7fffffffffffffffULL&hx; *eptr = 0; - if(ix>=0x7ff0000000000000ULL||((ix|ixl)==0)) return x; /* 0,inf,nan */ + if(ix>=0x7ff0000000000000ULL||ix==0) return x; /* 0,inf,nan */ if (ix<0x0010000000000000ULL) { /* subnormal */ x *= two107; - GET_LDOUBLE_MSW64(hx,x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); ix = hx&0x7fffffffffffffffULL; *eptr = -107; } @@ -54,7 +59,7 @@ long double __frexpl(long double x, int *eptr) if (ixl != 0ULL) { explo = (ixl>>52) - (ix>>52) + 0x3fe; if ((ixl&0x7ff0000000000000ULL) == 0LL) { - /* the lower double is a denomal so we need to correct its + /* the lower double is a denormal so we need to correct its mantissa and perhaps its exponent. */ int cnt; @@ -73,7 +78,9 @@ long double __frexpl(long double x, int *eptr) lx = 0ULL; hx = (hx&0x800fffffffffffffULL) | 0x3fe0000000000000ULL; - SET_LDOUBLE_WORDS64(x,hx,lx); + INSERT_WORDS64 (xhi, hx); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); return x; } #ifdef IS_IN_libm diff --git a/sysdeps/ieee754/ldbl-128ibm/s_isinf_nsl.c b/sysdeps/ieee754/ldbl-128ibm/s_isinf_nsl.c index c8dd9ff98a..54e72c9166 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_isinf_nsl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_isinf_nsl.c @@ -1,6 +1,7 @@ /* * __isinf_nsl(x) returns != 0 if x is ±inf, else 0; * no branching! + * slightly dodgy in relying on signed shift right copying sign bit */ #include <math.h> @@ -9,8 +10,14 @@ int __isinf_nsl (long double x) { - int64_t hx,lx; - GET_LDOUBLE_WORDS64(hx,lx,x); - return !((lx & 0x7fffffffffffffffLL) - | ((hx & 0x7fffffffffffffffLL) ^ 0x7ff0000000000000LL)); + double xhi; + int64_t hx, mask; + + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); + + mask = (hx & 0x7fffffffffffffffLL) ^ 0x7ff0000000000000LL; + mask |= -mask; + mask >>= 63; + return ~mask; } diff --git a/sysdeps/ieee754/ldbl-128ibm/s_isinfl.c b/sysdeps/ieee754/ldbl-128ibm/s_isinfl.c index 5f5b0144b2..6a728221fc 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_isinfl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_isinfl.c @@ -11,6 +11,7 @@ static char rcsid[] = "$NetBSD: $"; /* * isinfl(x) returns 1 if x is inf, -1 if x is -inf, else 0; * no branching! + * slightly dodgy in relying on signed shift right copying sign bit */ #include <math.h> @@ -20,12 +21,16 @@ static char rcsid[] = "$NetBSD: $"; int ___isinfl (long double x) { - int64_t hx,lx; - GET_LDOUBLE_WORDS64(hx,lx,x); - lx = (lx & 0x7fffffffffffffffLL); - lx |= (hx & 0x7fffffffffffffffLL) ^ 0x7ff0000000000000LL; - lx |= -lx; - return ~(lx >> 63) & (hx >> 62); + double xhi; + int64_t hx, mask; + + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); + + mask = (hx & 0x7fffffffffffffffLL) ^ 0x7ff0000000000000LL; + mask |= -mask; + mask >>= 63; + return ~mask & (hx >> 62); } hidden_ver (___isinfl, __isinfl) #ifndef IS_IN_libm diff --git a/sysdeps/ieee754/ldbl-128ibm/s_log1pl.c b/sysdeps/ieee754/ldbl-128ibm/s_log1pl.c index 77c4fdea84..a346383052 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_log1pl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_log1pl.c @@ -126,19 +126,18 @@ long double __log1pl (long double xm1) { long double x, y, z, r, s; - ieee854_long_double_shape_type u; - int32_t hx; + double xhi; + int32_t hx, lx; int e; /* Test for NaN or infinity input. */ - u.value = xm1; - hx = u.parts32.w0; + xhi = ldbl_high (xm1); + EXTRACT_WORDS (hx, lx, xhi); if (hx >= 0x7ff00000) return xm1; /* log1p(+- 0) = +- 0. */ - if (((hx & 0x7fffffff) == 0) - && (u.parts32.w1 | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3) == 0) + if (((hx & 0x7fffffff) | lx) == 0) return xm1; x = xm1 + 1.0L; diff --git a/sysdeps/ieee754/ldbl-128ibm/s_modfl.c b/sysdeps/ieee754/ldbl-128ibm/s_modfl.c index 39de9d4bfb..ed03ce236c 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_modfl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_modfl.c @@ -37,43 +37,54 @@ long double __modfl(long double x, long double *iptr) { int64_t i0,i1,j0; u_int64_t i; - GET_LDOUBLE_WORDS64(i0,i1,x); + double xhi, xlo; + + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (i0, xhi); + EXTRACT_WORDS64 (i1, xlo); i1 &= 0x000fffffffffffffLL; j0 = ((i0>>52)&0x7ff)-0x3ff; /* exponent of x */ if(j0<52) { /* integer part in high x */ if(j0<0) { /* |x|<1 */ /* *iptr = +-0 */ - SET_LDOUBLE_WORDS64(*iptr,i0&0x8000000000000000ULL,0); + INSERT_WORDS64 (xhi, i0&0x8000000000000000ULL); + *iptr = xhi; return x; } else { i = (0x000fffffffffffffLL)>>j0; if(((i0&i)|(i1&0x7fffffffffffffffLL))==0) { /* x is integral */ *iptr = x; /* return +-0 */ - SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0); + INSERT_WORDS64 (xhi, i0&0x8000000000000000ULL); + x = xhi; return x; } else { - SET_LDOUBLE_WORDS64(*iptr,i0&(~i),0); + INSERT_WORDS64 (xhi, i0&(~i)); + *iptr = xhi; return x - *iptr; } } } else if (j0>103) { /* no fraction part */ *iptr = x*one; /* We must handle NaNs separately. */ - if (j0 == 0x400 && ((i0 & 0x000fffffffffffffLL) | i1)) + if ((i0 & 0x7fffffffffffffffLL) > 0x7ff0000000000000LL) return x*one; /* return +-0 */ - SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0); + INSERT_WORDS64 (xhi, i0&0x8000000000000000ULL); + x = xhi; return x; } else { /* fraction part in low x */ i = -1ULL>>(j0-52); if((i1&i)==0) { /* x is integral */ *iptr = x; /* return +-0 */ - SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0); + INSERT_WORDS64 (xhi, i0&0x8000000000000000ULL); + x = xhi; return x; } else { - SET_LDOUBLE_WORDS64(*iptr,i0,i1&(~i)); + INSERT_WORDS64 (xhi, i0); + INSERT_WORDS64 (xlo, i1&(~i)); + *iptr = ldbl_pack (xhi, xlo); return x - *iptr; } } diff --git a/sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c b/sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c index 7e581274a3..c050944c0c 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c @@ -30,27 +30,28 @@ static char rcsid[] = "$NetBSD: $"; long double __nextafterl(long double x, long double y) { - int64_t hx,hy,ihx,ihy,ilx; - u_int64_t lx; - u_int64_t ly __attribute__ ((unused)); + int64_t hx,hy,ihx,ihy; + uint64_t lx; + double xhi, xlo, yhi; - GET_LDOUBLE_WORDS64(hx,lx,x); - GET_LDOUBLE_WORDS64(hy,ly,y); + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); + yhi = ldbl_high (y); + EXTRACT_WORDS64 (hy, yhi); ihx = hx&0x7fffffffffffffffLL; /* |hx| */ - ilx = lx&0x7fffffffffffffffLL; /* |lx| */ ihy = hy&0x7fffffffffffffffLL; /* |hy| */ - if((((ihx&0x7ff0000000000000LL)==0x7ff0000000000000LL)&& - ((ihx&0x000fffffffffffffLL)!=0)) || /* x is nan */ - (((ihy&0x7ff0000000000000LL)==0x7ff0000000000000LL)&& - ((ihy&0x000fffffffffffffLL)!=0))) /* y is nan */ + if((ihx>0x7ff0000000000000LL) || /* x is nan */ + (ihy>0x7ff0000000000000LL)) /* y is nan */ return x+y; /* signal the nan */ if(x==y) return y; /* x=y, return y */ - if(ihx == 0 && ilx == 0) { /* x == 0 */ - long double u; + if(ihx == 0) { /* x == 0 */ + long double u; /* return +-minsubnormal */ hy = (hy & 0x8000000000000000ULL) | 1; - SET_LDOUBLE_WORDS64(x,hy,0ULL);/* return +-minsubnormal */ + INSERT_WORDS64 (yhi, hy); + x = yhi; u = math_opt_barrier (x); u = u * u; math_force_eval (u); /* raise underflow flag */ @@ -59,10 +60,16 @@ long double __nextafterl(long double x, long double y) long double u; if(x > y) { /* x > y, x -= ulp */ + /* This isn't the largest magnitude correctly rounded + long double as you can see from the lowest mantissa + bit being zero. It is however the largest magnitude + long double with a 106 bit mantissa, and nextafterl + is insane with variable precision. So to make + nextafterl sane we assume 106 bit precision. */ if((hx==0xffefffffffffffffLL)&&(lx==0xfc8ffffffffffffeLL)) return x+x; /* overflow, return -inf */ if (hx >= 0x7ff0000000000000LL) { - SET_LDOUBLE_WORDS64(u,0x7fefffffffffffffLL,0x7c8ffffffffffffeLL); + u = 0x1.fffffffffffff7ffffffffffff8p+1023L; return u; } if(ihx <= 0x0360000000000000LL) { /* x <= LDBL_MIN */ @@ -77,16 +84,19 @@ long double __nextafterl(long double x, long double y) return x; } if (ihx < 0x06a0000000000000LL) { /* ulp will denormal */ - SET_LDOUBLE_WORDS64(u,(hx&0x7ff0000000000000LL),0ULL); + INSERT_WORDS64 (yhi, hx & (0x7ffLL<<52)); + u = yhi; u *= 0x1.0000000000000p-105L; - } else - SET_LDOUBLE_WORDS64(u,(hx&0x7ff0000000000000LL)-0x0690000000000000LL,0ULL); + } else { + INSERT_WORDS64 (yhi, (hx & (0x7ffLL<<52))-(0x069LL<<52)); + u = yhi; + } return x - u; } else { /* x < y, x += ulp */ if((hx==0x7fefffffffffffffLL)&&(lx==0x7c8ffffffffffffeLL)) return x+x; /* overflow, return +inf */ - if ((u_int64_t) hx >= 0xfff0000000000000ULL) { - SET_LDOUBLE_WORDS64(u,0xffefffffffffffffLL,0xfc8ffffffffffffeLL); + if ((uint64_t) hx >= 0xfff0000000000000ULL) { + u = -0x1.fffffffffffff7ffffffffffff8p+1023L; return u; } if(ihx <= 0x0360000000000000LL) { /* x <= LDBL_MIN */ @@ -103,10 +113,13 @@ long double __nextafterl(long double x, long double y) return x; } if (ihx < 0x06a0000000000000LL) { /* ulp will denormal */ - SET_LDOUBLE_WORDS64(u,(hx&0x7ff0000000000000LL),0ULL); + INSERT_WORDS64 (yhi, hx & (0x7ffLL<<52)); + u = yhi; u *= 0x1.0000000000000p-105L; - } else - SET_LDOUBLE_WORDS64(u,(hx&0x7ff0000000000000LL)-0x0690000000000000LL,0ULL); + } else { + INSERT_WORDS64 (yhi, (hx & (0x7ffLL<<52))-(0x069LL<<52)); + u = yhi; + } return x + u; } } diff --git a/sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c b/sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c index 7e288a4368..b40cf167f3 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c @@ -34,23 +34,22 @@ double __nexttoward(double x, long double y) { int32_t hx,ix; int64_t hy,iy; - u_int32_t lx; - u_int64_t ly,uly; + uint32_t lx; + double yhi; EXTRACT_WORDS(hx,lx,x); - GET_LDOUBLE_WORDS64(hy,ly,y); + yhi = ldbl_high (y); + EXTRACT_WORDS64(hy,yhi); ix = hx&0x7fffffff; /* |x| */ iy = hy&0x7fffffffffffffffLL; /* |y| */ - uly = ly&0x7fffffffffffffffLL; /* |y| */ if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) || /* x is nan */ - ((iy>=0x7ff0000000000000LL)&&((iy-0x7ff0000000000000LL)|uly)!=0)) - /* y is nan */ + iy>0x7ff0000000000000LL) /* y is nan */ return x+y; if((long double) x==y) return y; /* x=y, return y */ if((ix|lx)==0) { /* x == 0 */ double u; - INSERT_WORDS(x,(u_int32_t)((hy>>32)&0x80000000),1);/* return +-minsub */ + INSERT_WORDS(x,(uint32_t)((hy>>32)&0x80000000),1);/* return +-minsub */ u = math_opt_barrier (x); u = u * u; math_force_eval (u); /* raise underflow flag */ diff --git a/sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c b/sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c index b387a91192..19522f4762 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c @@ -27,16 +27,16 @@ float __nexttowardf(float x, long double y) { int32_t hx,ix; int64_t hy,iy; - u_int64_t ly, uly; + double yhi; GET_FLOAT_WORD(hx,x); - GET_LDOUBLE_WORDS64(hy,ly,y); + yhi = ldbl_high (y); + EXTRACT_WORDS64 (hy, yhi); ix = hx&0x7fffffff; /* |x| */ iy = hy&0x7fffffffffffffffLL; /* |y| */ - uly = ly&0x7fffffffffffffffLL; /* |y| */ if((ix>0x7f800000) || /* x is nan */ - ((iy>=0x7ff0000000000000LL)&&((iy-0x7ff0000000000000LL)|uly)!=0)) + (iy>0x7ff0000000000000LL)) /* y is nan */ return x+y; if((long double) x==y) return y; /* x=y, return y */ diff --git a/sysdeps/ieee754/ldbl-128ibm/s_remquol.c b/sysdeps/ieee754/ldbl-128ibm/s_remquol.c index f4777a0e1a..195e108ca9 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_remquol.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_remquol.c @@ -33,20 +33,24 @@ __remquol (long double x, long double y, int *quo) int64_t hx,hy; u_int64_t sx,lx,ly,qs; int cquo; - - GET_LDOUBLE_WORDS64 (hx, lx, x); - GET_LDOUBLE_WORDS64 (hy, ly, y); + double xhi, xlo, yhi, ylo; + + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); + ldbl_unpack (y, &yhi, &ylo); + EXTRACT_WORDS64 (hy, yhi); + EXTRACT_WORDS64 (ly, ylo); sx = hx & 0x8000000000000000ULL; qs = sx ^ (hy & 0x8000000000000000ULL); hy &= 0x7fffffffffffffffLL; hx &= 0x7fffffffffffffffLL; /* Purge off exception values. */ - if ((hy | (ly & 0x7fffffffffffffff)) == 0) + if (hy == 0) return (x * y) / (x * y); /* y = 0 */ if ((hx >= 0x7ff0000000000000LL) /* x not finite */ - || ((hy >= 0x7ff0000000000000LL) /* y is NaN */ - && (((hy - 0x7ff0000000000000LL) | ly) != 0))) + || (hy > 0x7ff0000000000000LL)) /* y is NaN */ return (x * y) / (x * y); if (hy <= 0x7fbfffffffffffffLL) diff --git a/sysdeps/ieee754/ldbl-128ibm/s_scalblnl.c b/sysdeps/ieee754/ldbl-128ibm/s_scalblnl.c index d752568885..03d4597271 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_scalblnl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_scalblnl.c @@ -41,11 +41,15 @@ long double __scalblnl (long double x, long int n) { int64_t k,l,hx,lx; union { int64_t i; double d; } u; - GET_LDOUBLE_WORDS64(hx,lx,x); + double xhi, xlo; + + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); k = (hx>>52)&0x7ff; /* extract exponent */ l = (lx>>52)&0x7ff; if (k==0) { /* 0 or subnormal x */ - if (((hx|lx)&0x7fffffffffffffffULL)==0) return x; /* +-0 */ + if ((hx&0x7fffffffffffffffULL)==0) return x; /* +-0 */ u.i = hx; u.d *= two54; hx = u.i; @@ -61,7 +65,9 @@ long double __scalblnl (long double x, long int n) if (k > 0) { /* normal result */ hx = (hx&0x800fffffffffffffULL)|(k<<52); if ((lx & 0x7fffffffffffffffULL) == 0) { /* low part +-0 */ - SET_LDOUBLE_WORDS64(x,hx,lx); + INSERT_WORDS64 (xhi, hx); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); return x; } if (l == 0) { /* low part subnormal */ @@ -81,14 +87,19 @@ long double __scalblnl (long double x, long int n) u.d *= twom54; lx = u.i; } - SET_LDOUBLE_WORDS64(x,hx,lx); + INSERT_WORDS64 (xhi, hx); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); return x; } if (k <= -54) return tiny*__copysignl(tiny,x); /*underflow*/ k += 54; /* subnormal result */ lx &= 0x8000000000000000ULL; - SET_LDOUBLE_WORDS64(x,(hx&0x800fffffffffffffULL)|(k<<52),lx); + hx &= 0x800fffffffffffffULL; + INSERT_WORDS64 (xhi, hx|(k<<52)); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); return x*twolm54; } long_double_symbol (libm, __scalblnl, scalblnl); diff --git a/sysdeps/ieee754/ldbl-128ibm/s_scalbnl.c b/sysdeps/ieee754/ldbl-128ibm/s_scalbnl.c index bcdb23bdaa..161172db6e 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_scalbnl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_scalbnl.c @@ -41,11 +41,15 @@ long double __scalbnl (long double x, int n) { int64_t k,l,hx,lx; union { int64_t i; double d; } u; - GET_LDOUBLE_WORDS64(hx,lx,x); + double xhi, xlo; + + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); k = (hx>>52)&0x7ff; /* extract exponent */ l = (lx>>52)&0x7ff; if (k==0) { /* 0 or subnormal x */ - if (((hx|lx)&0x7fffffffffffffffULL)==0) return x; /* +-0 */ + if ((hx&0x7fffffffffffffffULL)==0) return x; /* +-0 */ u.i = hx; u.d *= two54; hx = u.i; @@ -61,7 +65,9 @@ long double __scalbnl (long double x, int n) if (k > 0) { /* normal result */ hx = (hx&0x800fffffffffffffULL)|(k<<52); if ((lx & 0x7fffffffffffffffULL) == 0) { /* low part +-0 */ - SET_LDOUBLE_WORDS64(x,hx,lx); + INSERT_WORDS64 (xhi, hx); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); return x; } if (l == 0) { /* low part subnormal */ @@ -81,14 +87,19 @@ long double __scalbnl (long double x, int n) u.d *= twom54; lx = u.i; } - SET_LDOUBLE_WORDS64(x,hx,lx); + INSERT_WORDS64 (xhi, hx); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); return x; } if (k <= -54) return tiny*__copysignl(tiny,x); /*underflow*/ k += 54; /* subnormal result */ lx &= 0x8000000000000000ULL; - SET_LDOUBLE_WORDS64(x,(hx&0x800fffffffffffffULL)|(k<<52),lx); + hx &= 0x800fffffffffffffULL; + INSERT_WORDS64 (xhi, hx|(k<<52)); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); return x*twolm54; } #ifdef IS_IN_libm diff --git a/sysdeps/ieee754/ldbl-128ibm/s_tanhl.c b/sysdeps/ieee754/ldbl-128ibm/s_tanhl.c index 138b63cd1a..c63e25345d 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_tanhl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_tanhl.c @@ -47,10 +47,12 @@ static const long double one=1.0L, two=2.0L, tiny = 1.0e-300L; long double __tanhl(long double x) { long double t,z; - int64_t jx,ix,lx; + int64_t jx,ix; + double xhi; /* High word of |x|. */ - GET_LDOUBLE_WORDS64(jx,lx,x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (jx, xhi); ix = jx&0x7fffffffffffffffLL; /* x is INF or NaN */ @@ -61,7 +63,7 @@ long double __tanhl(long double x) /* |x| < 22 */ if (ix < 0x4036000000000000LL) { /* |x|<22 */ - if ((ix | (lx&0x7fffffffffffffffLL)) == 0) + if (ix == 0) return x; /* x == +-0 */ if (ix<0x3c60000000000000LL) /* |x|<2**-57 */ return x*(one+x); /* tanh(small) = small */ diff --git a/sysdeps/powerpc/fpu/libm-test-ulps b/sysdeps/powerpc/fpu/libm-test-ulps index 74365f02fa..37b2ca192d 100644 --- a/sysdeps/powerpc/fpu/libm-test-ulps +++ b/sysdeps/powerpc/fpu/libm-test-ulps @@ -6640,6 +6640,9 @@ float: 1 ifloat: 1 ildouble: 2 ldouble: 2 +Test "tan_towardzero (2)": +ildouble: 1 +ldouble: 1 Test "tan_towardzero (3)": float: 1 ifloat: 1 |