* [PATCH] Fix a couple of ldbl-128 bugs
@ 2002-07-03 3:52 Jakub Jelinek
2002-07-10 20:13 ` Ulrich Drepper
0 siblings, 1 reply; 2+ messages in thread
From: Jakub Jelinek @ 2002-07-03 3:52 UTC (permalink / raw)
To: Ulrich Drepper, Stephen L Moshier; +Cc: Glibc hackers
Hi!
While trying to test glibc sizeof(long double) == 8, sizeof(long double) == 16
patch (depending on compiler switches) I've found the following bugs.
printf was broken with some numbers, because
if _FPIO_CONST_SHIFT was non-zero (only IEEE quad long double on
32-bit platforms), scaleexpo would be increased too much (64bits more) and
thus if (exponent >= scaleexpo + powers->p_expo - 1)
line would never trigger any more. This had the result that scale could be
way more than 10 times smaller than frac, resulting in the first digit
of the number (char)('0' + some huge number), e.g. printing
o01757308.0000
instead of
43679901757308.0000
#ifdef NO_LONG_DOUBLE aliases make no sense in ldbl-128 directory, as using
that directory implies long double is supported.
There were several bugs in __mpn_extract_long_double for
32-bit arches, remquol was broken because qs was int and doing
qs = sx ^ (hy & 0x8000000000000000ULL);
implied qs was always 0 due to limited type range.
_FP_FRAC_CLZ_4 had a pasto which caused bad handling of quad denormals
on 32-bit arches.
And last this patch fixes a bunch of test-ldouble failures where
either exceptions were thrown where they shouldn't or
vice versa (I've tried to follow what dbl-64 and flt-32 implementations do
for the special cases).
2002-07-03 Jakub Jelinek <jakub@redhat.com>
* stdio-common/printf_fp.c (__printf_fp.c): If _FPIO_CONST_SHIFT is
non-zero, adjust exponent.
* sysdeps/ieee754/ldbl-128/s_erfl.c (__erfl, erfl, __erfcl, erfcl):
Remove NO_LONG_DOUBLE aliases.
* sysdeps/ieee754/ldbl-128/s_expm1l.c (__expm1l, expm1l): Likewise.
* sysdeps/ieee754/ldbl-128/s_log1pl.c (__log1pl, log1pl): Likewise.
(__log1pl): Raise divide by zero and invalid exceptions when needed.
* sysdeps/ieee754/ldbl-128/e_powl.c (__ieee754_powl): Special case
1**y and -1**+-Inf.
* sysdeps/ieee754/ldbl-128/ldbl2mpn.c (__mpn_extract_long_double):
Fix BITS_PER_MP_LIMB 32 extraction.
* sysdeps/ieee754/ldbl-128/e_log2l.c (__ieee754_log2l): Don't raise
exceptions for qNaNs.
* sysdeps/ieee754/ldbl-128/e_log10l.c (__ieee754_log10l): Likewise.
* sysdeps/ieee754/ldbl-128/e_lgammal_r.c (__ieee754_lgamma_r):
Raise exceptions when needed. Don't recurse unnecessarily.
Special case 1.0L and 2.0L arguments to avoid -0.0L as result.
* sysdeps/ieee754/ldbl-128/e_j0l.c (__ieee754_y0l): Don't raise
exceptions for qNaNs.
* sysdeps/ieee754/ldbl-128/s_remquol.c (__remquol): Make qs 64-bit
to fix *quo return value sign.
* sysdeps/ieee754/ldbl-128/e_gammal_r.c (__ieee754_gamma_r): Special
case -Inf argument.
* soft-fp/op-4.h (_FP_FRAC_CLZ_4): Fix a pasto.
--- libc/stdio-common/printf_fp.c.jj Tue Apr 30 12:53:29 2002
+++ libc/stdio-common/printf_fp.c Tue Jul 2 15:35:23 2002
@@ -494,6 +494,9 @@ __printf_fp (FILE *fp,
&__tens[powers->arrayoff],
tmpsize * sizeof (mp_limb_t));
MPN_ZERO (tmp, _FPIO_CONST_SHIFT);
+ /* Adjust exponent, as scaleexpo will be this much
+ bigger too. */
+ exponent += _FPIO_CONST_SHIFT * BITS_PER_MP_LIMB;
}
else
#endif
--- libc/sysdeps/ieee754/ldbl-128/s_erfl.c.jj Thu Nov 29 13:38:01 2001
+++ libc/sysdeps/ieee754/ldbl-128/s_erfl.c Mon Jul 1 12:56:48 2002
@@ -790,10 +790,6 @@ __erfl (x)
}
weak_alias (__erfl, erfl)
-#ifdef NO_LONG_DOUBLE
-strong_alias (__erf, __erfl)
-weak_alias (__erf, erfl)
-#endif
#ifdef __STDC__
long double
__erfcl (long double x)
@@ -935,7 +931,3 @@ weak_alias (__erf, erfl)
}
weak_alias (__erfcl, erfcl)
-#ifdef NO_LONG_DOUBLE
-strong_alias (__erfc, __erfcl)
-weak_alias (__erfc, erfcl)
-#endif
--- libc/sysdeps/ieee754/ldbl-128/e_powl.c.jj Sun Oct 14 23:38:42 2001
+++ libc/sysdeps/ieee754/ldbl-128/e_powl.c Tue Jul 2 18:24:04 2002
@@ -156,6 +156,13 @@ __ieee754_powl (x, y)
if ((iy | q.parts32.w1 | q.parts32.w2 | q.parts32.w3) == 0)
return one;
+ /* 1.0**y = 1; -1.0**+-Inf = 1 */
+ if (x == one)
+ return one;
+ if (x == -1.0L && iy == 0x7fff0000
+ && (q.parts32.w1 | q.parts32.w2 | q.parts32.w3) == 0)
+ return one;
+
/* +-NaN return x+y */
if ((ix > 0x7fff0000)
|| ((ix == 0x7fff0000)
--- libc/sysdeps/ieee754/ldbl-128/ldbl2mpn.c.jj Thu Aug 23 18:49:59 2001
+++ libc/sysdeps/ieee754/ldbl-128/ldbl2mpn.c Mon Jul 1 20:23:27 2002
@@ -102,7 +102,7 @@ __mpn_extract_long_double (mp_ptr res_pt
#else
int j, k, l;
- for (j = N - 1; j > 0; j++)
+ for (j = N - 1; j > 0; j--)
if (res_ptr[j] != 0)
break;
@@ -112,20 +112,22 @@ __mpn_extract_long_double (mp_ptr res_pt
if (cnt < 0)
{
cnt += BITS_PER_MP_LIMB;
- l++;
+ l--;
}
if (!cnt)
for (k = N - 1; k >= l; k--)
res_ptr[k] = res_ptr[k-l];
else
- for (k = N - 1; k >= l; k--)
- res_ptr[k] = res_ptr[k-l] << cnt
- | res_ptr[k-l-1] >> (BITS_PER_MP_LIMB - cnt);
- res_ptr[k--] = res_ptr[0] << cnt;
+ {
+ for (k = N - 1; k > l; k--)
+ res_ptr[k] = res_ptr[k-l] << cnt
+ | res_ptr[k-l-1] >> (BITS_PER_MP_LIMB - cnt);
+ res_ptr[k--] = res_ptr[0] << cnt;
+ }
for (; k >= 0; k--)
res_ptr[k] = 0;
- *expt = LDBL_MIN_EXP - 1 - 3 * BITS_PER_MP_LIMB - cnt;
+ *expt = LDBL_MIN_EXP - 1 - l * BITS_PER_MP_LIMB - cnt;
#endif
}
}
--- libc/sysdeps/ieee754/ldbl-128/e_log2l.c.jj Sat Nov 10 11:38:27 2001
+++ libc/sysdeps/ieee754/ldbl-128/e_log2l.c Tue Jul 2 16:44:22 2002
@@ -164,16 +164,15 @@ __ieee754_log2l (x)
long double z;
long double y;
int e;
+ int64_t hx, lx;
/* Test for domain */
- if (x <= 0.0L)
- {
- if (x == 0.0L)
- return (-1.0L / (x - x));
- else
- return (x - x) / (x - x);
- }
- if (!__finitel (x))
+ GET_LDOUBLE_WORDS64 (hx, lx, x);
+ if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
+ return (-1.0L / (x - x));
+ if (hx < 0)
+ return (x - x) / (x - x);
+ if (hx >= 0x7fff000000000000LL)
return (x + x);
/* separate mantissa from exponent */
--- libc/sysdeps/ieee754/ldbl-128/e_log10l.c.jj Thu Nov 29 13:38:01 2001
+++ libc/sysdeps/ieee754/ldbl-128/e_log10l.c Tue Jul 2 16:44:54 2002
@@ -170,16 +170,15 @@ __ieee754_log10l (x)
long double z;
long double y;
int e;
+ int64_t hx, lx;
/* Test for domain */
- if (x <= 0.0L)
- {
- if (x == 0.0L)
- return (-1.0L / (x - x));
- else
- return (x - x) / (x - x);
- }
- if (!__finitel (x))
+ GET_LDOUBLE_WORDS64 (hx, lx, x);
+ if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
+ return (-1.0L / (x - x));
+ if (hx < 0)
+ return (x - x) / (x - x);
+ if (hx >= 0x7fff000000000000LL)
return (x + x);
/* separate mantissa from exponent */
--- libc/sysdeps/ieee754/ldbl-128/s_expm1l.c.jj Wed Nov 21 13:32:44 2001
+++ libc/sysdeps/ieee754/ldbl-128/s_expm1l.c Mon Jul 1 12:56:48 2002
@@ -144,6 +144,3 @@ __expm1l (long double x)
}
weak_alias (__expm1l, expm1l)
-#ifdef NO_LONG_DOUBLE
-strong_alias (__expm1, __expm1l) weak_alias (__expm1, expm1l)
-#endif
--- libc/sysdeps/ieee754/ldbl-128/s_log1pl.c.jj Wed Nov 21 13:32:44 2001
+++ libc/sysdeps/ieee754/ldbl-128/s_log1pl.c Tue Jul 2 17:37:53 2002
@@ -117,17 +117,18 @@ __log1pl (long double xm1)
{
long double x, y, z, r, s;
ieee854_long_double_shape_type u;
- int32_t ix;
+ int32_t hx;
int e;
/* Test for NaN or infinity input. */
u.value = xm1;
- ix = u.parts32.w0 & 0x7fffffff;
- if (ix >= 0x7fff0000)
+ hx = u.parts32.w0;
+ if (hx >= 0x7fff0000)
return xm1;
/* log1p(+- 0) = +- 0. */
- if ((ix == 0) && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
+ if (((hx & 0x7fffffff) == 0)
+ && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
return xm1;
x = xm1 + 1.0L;
@@ -136,9 +137,9 @@ __log1pl (long double xm1)
if (x <= 0.0L)
{
if (x == 0.0L)
- return (-1.0L / zero);
+ return (-1.0L / (x - x));
else
- return (zero / zero);
+ return (zero / (x - x));
}
/* Separate mantissa from exponent. */
@@ -238,7 +239,3 @@ __log1pl (long double xm1)
}
weak_alias (__log1pl, log1pl)
-#ifdef NO_LONG_DOUBLE
-strong_alias (__log1p, __log1pl)
-weak_alias (__log1p, log1pl)
-#endif
--- libc/sysdeps/ieee754/ldbl-128/e_lgammal_r.c.jj Thu Nov 29 13:38:01 2001
+++ libc/sysdeps/ieee754/ldbl-128/e_lgammal_r.c Tue Jul 2 22:29:26 2002
@@ -761,10 +761,9 @@ __ieee754_lgammal_r (x, signgamp)
if (x < 0.0L)
{
q = -x;
- w = __ieee754_lgammal_r (q, &i);
p = __floorl (q);
if (p == q)
- return (one / zero);
+ return (one / (p - p));
i = p;
if ((i & 1) == 0)
*signgamp = -1;
@@ -779,6 +778,7 @@ __ieee754_lgammal_r (x, signgamp)
z = q * __sinl (PIL * z);
if (z == 0.0L)
return (*signgamp * huge * huge);
+ w = __ieee754_lgammal_r (q, &i);
z = __logl (PIL / z) - w;
return (z);
}
@@ -859,6 +859,8 @@ __ieee754_lgammal_r (x, signgamp)
z = x - 1.0L;
p = z * neval (z, RNr9, NRNr9) / deval (z, RDr9, NRDr9);
}
+ else if (x == 1.0L)
+ p = 0.0L;
else if (x <= 1.125L)
{
z = x - 1.0L;
@@ -900,6 +902,8 @@ __ieee754_lgammal_r (x, signgamp)
p += lgam1r75b;
p += lgam1r75a;
}
+ else if (x == 2.0L)
+ p = 0.0L;
else if (x < 2.375L)
{
z = x - 2.0L;
--- libc/sysdeps/ieee754/ldbl-128/e_j0l.c.jj Thu Nov 29 13:38:01 2001
+++ libc/sysdeps/ieee754/ldbl-128/e_j0l.c Tue Jul 2 19:33:10 2002
@@ -803,12 +803,6 @@ long double
{
long double xx, xinv, z, p, q, c, s, cc, ss;
- if (x <= 0.0L)
- {
- if (x < 0.0L)
- return (zero / zero);
- return 1.0L / zero;
- }
if (! finitel (x))
{
if (x != x)
@@ -816,6 +810,12 @@ long double
else
return 0.0L;
}
+ if (x <= 0.0L)
+ {
+ if (x < 0.0L)
+ return (zero / zero);
+ return 1.0L / zero;
+ }
xx = fabsl (x);
if (xx <= 2.0L)
{
--- libc/sysdeps/ieee754/ldbl-128/s_remquol.c.jj Thu Aug 23 18:49:59 2001
+++ libc/sysdeps/ieee754/ldbl-128/s_remquol.c Tue Jul 2 20:15:00 2002
@@ -1,5 +1,5 @@
/* Compute remainder and a congruent to the quotient.
- Copyright (C) 1997, 1999 Free Software Foundation, Inc.
+ Copyright (C) 1997, 1999, 2002 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
Jakub Jelinek <jj@ultra.linux.cz>, 1999.
@@ -31,8 +31,8 @@ long double
__remquol (long double x, long double y, int *quo)
{
int64_t hx,hy;
- u_int64_t sx,lx,ly;
- int cquo,qs;
+ u_int64_t sx,lx,ly,qs;
+ int cquo;
GET_LDOUBLE_WORDS64 (hx, lx, x);
GET_LDOUBLE_WORDS64 (hy, ly, y);
--- libc/sysdeps/ieee754/ldbl-128/e_gammal_r.c.jj Thu Aug 23 18:49:54 2001
+++ libc/sysdeps/ieee754/ldbl-128/e_gammal_r.c Tue Jul 2 21:52:38 2002
@@ -1,5 +1,5 @@
/* Implementation of gamma function according to ISO C.
- Copyright (C) 1997, 1999 Free Software Foundation, Inc.
+ Copyright (C) 1997, 1999, 2002 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
Jakub Jelinek <jj@ultra.linux.cz, 1999.
@@ -46,6 +46,12 @@ __ieee754_gammal_r (long double x, int *
*signgamp = 0;
return (x - x) / (x - x);
}
+ if (hx == 0xffff000000000000ULL && lx == 0)
+ {
+ /* x == -Inf. According to ISO this is NaN. */
+ *signgamp = 0;
+ return x - x;
+ }
/* XXX FIXME. */
return __ieee754_expl (__ieee754_lgammal_r (x, signgamp));
--- libc/soft-fp/op-4.h.jj Tue Jul 2 09:46:19 2002
+++ libc/soft-fp/op-4.h Tue Jul 2 09:46:19 2002
@@ -167,7 +167,7 @@
} \
else if (X##_f[1]) \
{ \
- __FP_CLZ(R,X##_f[2]); \
+ __FP_CLZ(R,X##_f[1]); \
R += _FP_W_TYPE_SIZE*2; \
} \
else \
Jakub
^ permalink raw reply [flat|nested] 2+ messages in thread
* Re: [PATCH] Fix a couple of ldbl-128 bugs
2002-07-03 3:52 [PATCH] Fix a couple of ldbl-128 bugs Jakub Jelinek
@ 2002-07-10 20:13 ` Ulrich Drepper
0 siblings, 0 replies; 2+ messages in thread
From: Ulrich Drepper @ 2002-07-10 20:13 UTC (permalink / raw)
To: Jakub Jelinek; +Cc: Stephen L Moshier, Glibc hackers
[-- Attachment #1: Type: text/plain, Size: 545 bytes --]
On Wed, 2002-07-03 at 03:52, Jakub Jelinek wrote:
> While trying to test glibc sizeof(long double) == 8, sizeof(long double) == 16
> patch (depending on compiler switches) I've found the following bugs.
I have no way to test the patch but you certainly did and the patch
looks fine. I've applied it. Thanks,
--
---------------. ,-. 1325 Chesapeake Terrace
Ulrich Drepper \ ,-------------------' \ Sunnyvale, CA 94089 USA
Red Hat `--' drepper at redhat.com `------------------------
[-- Attachment #2: This is a digitally signed message part --]
[-- Type: application/pgp-signature, Size: 232 bytes --]
^ permalink raw reply [flat|nested] 2+ messages in thread
end of thread, other threads:[~2002-07-11 3:13 UTC | newest]
Thread overview: 2+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2002-07-03 3:52 [PATCH] Fix a couple of ldbl-128 bugs Jakub Jelinek
2002-07-10 20:13 ` Ulrich Drepper
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).