mirror of
https://sourceware.org/git/glibc.git
synced 2024-12-26 04:31:03 +00:00
cc4084017e
The ldbl-128ibm implementation of remainderl has logic resulting in incorrect tests for equality of the absolute values of the arguments in the case of zero low parts. If the low parts are both zero but with different signs, this can wrongly cause equal arguments to be treated as different, resulting in turn in incorrect signs of zero result in nondefault rounding modes arising from the subtractions done when the arguments are not equal. This patch fixes the logic to convert -0 low parts into +0 before the comparison (remquo already has separate logic to deal with signs of zero results, so doesn't need such a change). Tests are added for remainderl and remquol similar to that for fmodl, and based on a refactoring of it, since the bug depends on low parts which should not be relied upon in tests not setting the representation explicitly (although in fact the bug shows up in test-ldouble with current GCC). Tested for powerpc. [BZ #19677] * sysdeps/ieee754/ldbl-128ibm/e_remainderl.c (__ieee754_remainderl): Put zero low parts in canonical form. * sysdeps/ieee754/ldbl-128ibm/test-fmodrem-ldbl-128ibm.c: New file. Based on sysdeps/ieee754/ldbl-128ibm/test-fmodl-ldbl-128ibm.c. * sysdeps/ieee754/ldbl-128ibm/test-fmodl-ldbl-128ibm.c: Replace with wrapper round test-fmodrem-ldbl-128ibm.c. * sysdeps/ieee754/ldbl-128ibm/test-remainderl-ldbl-128ibm.c: New file. * sysdeps/ieee754/ldbl-128ibm/test-remquol-ldbl-128ibm.c: Likewise. * sysdeps/ieee754/ldbl-128ibm/Makefile (tests): Add test-remainderl-ldbl-128ibm and test-remquol-ldbl-128ibm.
82 lines
2.0 KiB
C
82 lines
2.0 KiB
C
/* e_fmodl.c -- long double version of e_fmod.c.
|
|
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
|
|
*/
|
|
/*
|
|
* ====================================================
|
|
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
|
*
|
|
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
|
* Permission to use, copy, modify, and distribute this
|
|
* software is freely granted, provided that this notice
|
|
* is preserved.
|
|
* ====================================================
|
|
*/
|
|
|
|
/* __ieee754_remainderl(x,p)
|
|
* Return :
|
|
* returns x REM p = x - [x/p]*p as if in infinite
|
|
* precise arithmetic, where [x/p] is the (infinite bit)
|
|
* integer nearest x/p (in half way case choose the even one).
|
|
* Method :
|
|
* Based on fmodl() return x-[x/p]chopped*p exactlp.
|
|
*/
|
|
|
|
#include <math.h>
|
|
#include <math_private.h>
|
|
|
|
static const long double zero = 0.0L;
|
|
|
|
|
|
long double
|
|
__ieee754_remainderl(long double x, long double p)
|
|
{
|
|
int64_t hx,hp;
|
|
u_int64_t sx,lx,lp;
|
|
long double p_half;
|
|
double xhi, xlo, phi, plo;
|
|
|
|
ldbl_unpack (x, &xhi, &xlo);
|
|
EXTRACT_WORDS64 (hx, xhi);
|
|
EXTRACT_WORDS64 (lx, xlo);
|
|
ldbl_unpack (p, &phi, &plo);
|
|
EXTRACT_WORDS64 (hp, phi);
|
|
EXTRACT_WORDS64 (lp, plo);
|
|
sx = hx&0x8000000000000000ULL;
|
|
lp ^= hp & 0x8000000000000000ULL;
|
|
hp &= 0x7fffffffffffffffLL;
|
|
lx ^= sx;
|
|
hx &= 0x7fffffffffffffffLL;
|
|
if (lp == 0x8000000000000000ULL)
|
|
lp = 0;
|
|
if (lx == 0x8000000000000000ULL)
|
|
lx = 0;
|
|
|
|
/* purge off exception values */
|
|
if(hp==0) return (x*p)/(x*p); /* p = 0 */
|
|
if((hx>=0x7ff0000000000000LL)|| /* x not finite */
|
|
(hp>0x7ff0000000000000LL)) /* p is NaN */
|
|
return (x*p)/(x*p);
|
|
|
|
|
|
if (hp<=0x7fdfffffffffffffLL) x = __ieee754_fmodl(x,p+p); /* now x < 2p */
|
|
if (((hx-hp)|(lx-lp))==0) return zero*x;
|
|
x = fabsl(x);
|
|
p = fabsl(p);
|
|
if (hp<0x0020000000000000LL) {
|
|
if(x+x>p) {
|
|
x-=p;
|
|
if(x+x>=p) x -= p;
|
|
}
|
|
} else {
|
|
p_half = 0.5L*p;
|
|
if(x>p_half) {
|
|
x-=p;
|
|
if(x>=p_half) x -= p;
|
|
}
|
|
}
|
|
if (sx)
|
|
x = -x;
|
|
return x;
|
|
}
|
|
strong_alias (__ieee754_remainderl, __remainderl_finite)
|