d9afe48d55
The dbl-64/wordsize-64 remquo implementation follows similar logic to various other implementations, but where that logic computes some absolute values, it wrongly uses a previously computed bit-pattern for the absolute value of the first argument, where actually it needs the absolute value of the first argument mod 8 times the second. This patch fixes it to compute the correct absolute value. The integer quotient result of remquo is only specified mod 8 (including its sign); architecture-specific versions may well vary in what results they give for higher bits of that result (and indeed bug 17569 gives an example correct result from __builtin_remquo giving 9 for that result, where the particular glibc implementation used in that bug report would give 1 after this fix). Thus, this patch adapts the tests of remquo to test that result only mod 8, to allow for such variation when tests with higher quotient are included. Tested for x86_64 and x86. [BZ #17569] * sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c (__remquo): Compute absolute value of x as modified by fmod, not original value of x. * math/libm-test.inc (RUN_TEST_ffI_f1): Rename to RUN_TEST_ffI_f1_mod8. Check extra return value mod 8. (RUN_TEST_LOOP_ffI_f1): Rename to RUN_TEST_LOOP_ffI_f1_mod8. Call RUN_TEST_ffI_f1_mod8. (remquo_test_data): Add more tests.
112 lines
2.4 KiB
C
112 lines
2.4 KiB
C
/* Compute remainder and a congruent to the quotient.
|
|
Copyright (C) 1997-2015 Free Software Foundation, Inc.
|
|
This file is part of the GNU C Library.
|
|
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
|
|
|
|
The GNU C Library is free software; you can redistribute it and/or
|
|
modify it under the terms of the GNU Lesser General Public
|
|
License as published by the Free Software Foundation; either
|
|
version 2.1 of the License, or (at your option) any later version.
|
|
|
|
The GNU C Library is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Lesser General Public License for more details.
|
|
|
|
You should have received a copy of the GNU Lesser General Public
|
|
License along with the GNU C Library; if not, see
|
|
<http://www.gnu.org/licenses/>. */
|
|
|
|
#include <math.h>
|
|
|
|
#include <math_private.h>
|
|
#include <stdint.h>
|
|
|
|
static const double zero = 0.0;
|
|
|
|
|
|
double
|
|
__remquo (double x, double y, int *quo)
|
|
{
|
|
int64_t hx, hy;
|
|
uint64_t sx, qs;
|
|
int cquo;
|
|
|
|
EXTRACT_WORDS64 (hx, x);
|
|
EXTRACT_WORDS64 (hy, y);
|
|
sx = hx & UINT64_C(0x8000000000000000);
|
|
qs = sx ^ (hy & UINT64_C(0x8000000000000000));
|
|
hy &= UINT64_C(0x7fffffffffffffff);
|
|
hx &= UINT64_C(0x7fffffffffffffff);
|
|
|
|
/* Purge off exception values. */
|
|
if (__glibc_unlikely (hy == 0))
|
|
return (x * y) / (x * y); /* y = 0 */
|
|
if (__builtin_expect (hx >= UINT64_C(0x7ff0000000000000) /* x not finite */
|
|
|| hy > UINT64_C(0x7ff0000000000000), 0))/* y is NaN */
|
|
return (x * y) / (x * y);
|
|
|
|
if (hy <= UINT64_C(0x7fbfffffffffffff))
|
|
x = __ieee754_fmod (x, 8 * y); /* now x < 8y */
|
|
|
|
if (__glibc_unlikely (hx == hy))
|
|
{
|
|
*quo = qs ? -1 : 1;
|
|
return zero * x;
|
|
}
|
|
|
|
x = fabs (x);
|
|
INSERT_WORDS64 (y, hy);
|
|
cquo = 0;
|
|
|
|
if (x >= 4 * y)
|
|
{
|
|
x -= 4 * y;
|
|
cquo += 4;
|
|
}
|
|
if (x >= 2 * y)
|
|
{
|
|
x -= 2 * y;
|
|
cquo += 2;
|
|
}
|
|
|
|
if (hy < UINT64_C(0x0020000000000000))
|
|
{
|
|
if (x + x > y)
|
|
{
|
|
x -= y;
|
|
++cquo;
|
|
if (x + x >= y)
|
|
{
|
|
x -= y;
|
|
++cquo;
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
double y_half = 0.5 * y;
|
|
if (x > y_half)
|
|
{
|
|
x -= y;
|
|
++cquo;
|
|
if (x >= y_half)
|
|
{
|
|
x -= y;
|
|
++cquo;
|
|
}
|
|
}
|
|
}
|
|
|
|
*quo = qs ? -cquo : cquo;
|
|
|
|
if (sx)
|
|
x = -x;
|
|
return x;
|
|
}
|
|
weak_alias (__remquo, remquo)
|
|
#ifdef NO_LONG_DOUBLE
|
|
strong_alias (__remquo, __remquol)
|
|
weak_alias (__remquo, remquol)
|
|
#endif
|