2011-09-08 02:10:26 +00:00
|
|
|
/* Pythagorean addition using doubles
|
2017-01-01 00:14:16 +00:00
|
|
|
Copyright (C) 2011-2017 Free Software Foundation, Inc.
|
2011-09-08 02:10:26 +00:00
|
|
|
This file is part of the GNU C Library
|
|
|
|
Contributed by Adhemerval Zanella <azanella@br.ibm.com>, 2011
|
|
|
|
|
|
|
|
The GNU C Library is free software; you can redistribute it and/or
|
|
|
|
modify it under the terms of the GNU Library General Public License as
|
|
|
|
published by the Free Software Foundation; either version 2 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
|
|
|
|
Library General Public License for more details.
|
|
|
|
|
|
|
|
You should have received a copy of the GNU Library General Public
|
2012-02-09 23:18:22 +00:00
|
|
|
License along with the GNU C Library; see the file COPYING.LIB. If
|
|
|
|
not, see <http://www.gnu.org/licenses/>. */
|
2011-09-08 02:10:26 +00:00
|
|
|
|
2012-03-09 19:29:16 +00:00
|
|
|
#include <math.h>
|
|
|
|
#include <math_private.h>
|
2013-05-01 15:46:34 +00:00
|
|
|
#include <stdint.h>
|
2011-09-08 02:10:26 +00:00
|
|
|
|
|
|
|
static const double two60 = 1.152921504606847e+18;
|
|
|
|
static const double two500 = 3.2733906078961419e+150;
|
|
|
|
static const double two600 = 4.149515568880993e+180;
|
|
|
|
static const double two1022 = 4.49423283715579e+307;
|
|
|
|
static const double twoM500 = 3.054936363499605e-151;
|
2011-12-06 10:10:06 +00:00
|
|
|
static const double twoM600 = 2.4099198651028841e-181;
|
2013-05-06 19:40:17 +00:00
|
|
|
static const double two60factor = 1.5592502418239997e+290;
|
2011-09-08 02:10:26 +00:00
|
|
|
static const double pdnum = 2.225073858507201e-308;
|
|
|
|
|
|
|
|
/* __ieee754_hypot(x,y)
|
|
|
|
*
|
|
|
|
* This a FP only version without any FP->INT conversion.
|
|
|
|
* It is similar to default C version, making appropriates
|
|
|
|
* overflow and underflows checks as well scaling when it
|
|
|
|
* is needed.
|
|
|
|
*/
|
|
|
|
|
|
|
|
#ifdef _ARCH_PWR7
|
|
|
|
/* POWER7 isinf and isnan optimization are fast. */
|
|
|
|
# define TEST_INF_NAN(x, y) \
|
Fix hypot sNaN handling (bug 20940).
TS 18661-1 generally defines libm functions taking sNaN arguments to
return qNaN and raise "invalid", even for the cases where a
corresponding qNaN argument would not result in a qNaN return. This
includes hypot with one argument being an infinity and the other being
an sNaN. This patch duly fixes hypot implementatations in glibc
(generic and powerpc) to ensure qNaN, computed by arithmetic on the
arguments, is returned in that case.
Various implementations do their checks for infinities and NaNs inline
by manipulating the representations of the arguments. For simplicity,
this patch just uses issignaling to check for sNaN arguments. This
could be inlined like the existing code (with due care about reversed
quiet NaN conventions, for implementations where that is relevant),
but given that all these checks are in cases where it's already known
at least one argument is not finite, which should be the uncommon
case, that doesn't seem worthwhile unless performance issues are
observed in practice.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #20940]
* sysdeps/ieee754/dbl-64/e_hypot.c (__ieee754_hypot): Do not
return Inf for arguments Inf and sNaN.
* sysdeps/ieee754/flt-32/e_hypotf.c (__ieee754_hypotf): Likewise.
* sysdeps/ieee754/ldbl-128/e_hypotl.c (__ieee754_hypotl):
Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_hypotl.c (__ieee754_hypotl):
Likewise.
* sysdeps/ieee754/ldbl-96/e_hypotl.c (__ieee754_hypotl): Likewise.
* sysdeps/powerpc/fpu/e_hypot.c (TEST_INF_NAN): Do not return Inf
for arguments Inf and sNaN. When returning a NaN, compute it by
arithmetic on the arguments.
* sysdeps/powerpc/fpu/e_hypotf.c (TEST_INF_NAN): Likewise.
* math/libm-test.inc (pow_test_data): Add tests of sNaN arguments.
2016-12-07 01:16:36 +00:00
|
|
|
if ((isinf(x) || isinf(y)) \
|
|
|
|
&& !issignaling (x) && !issignaling (y)) \
|
2011-09-08 02:10:26 +00:00
|
|
|
return INFINITY; \
|
|
|
|
if (isnan(x) || isnan(y)) \
|
Fix hypot sNaN handling (bug 20940).
TS 18661-1 generally defines libm functions taking sNaN arguments to
return qNaN and raise "invalid", even for the cases where a
corresponding qNaN argument would not result in a qNaN return. This
includes hypot with one argument being an infinity and the other being
an sNaN. This patch duly fixes hypot implementatations in glibc
(generic and powerpc) to ensure qNaN, computed by arithmetic on the
arguments, is returned in that case.
Various implementations do their checks for infinities and NaNs inline
by manipulating the representations of the arguments. For simplicity,
this patch just uses issignaling to check for sNaN arguments. This
could be inlined like the existing code (with due care about reversed
quiet NaN conventions, for implementations where that is relevant),
but given that all these checks are in cases where it's already known
at least one argument is not finite, which should be the uncommon
case, that doesn't seem worthwhile unless performance issues are
observed in practice.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #20940]
* sysdeps/ieee754/dbl-64/e_hypot.c (__ieee754_hypot): Do not
return Inf for arguments Inf and sNaN.
* sysdeps/ieee754/flt-32/e_hypotf.c (__ieee754_hypotf): Likewise.
* sysdeps/ieee754/ldbl-128/e_hypotl.c (__ieee754_hypotl):
Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_hypotl.c (__ieee754_hypotl):
Likewise.
* sysdeps/ieee754/ldbl-96/e_hypotl.c (__ieee754_hypotl): Likewise.
* sysdeps/powerpc/fpu/e_hypot.c (TEST_INF_NAN): Do not return Inf
for arguments Inf and sNaN. When returning a NaN, compute it by
arithmetic on the arguments.
* sysdeps/powerpc/fpu/e_hypotf.c (TEST_INF_NAN): Likewise.
* math/libm-test.inc (pow_test_data): Add tests of sNaN arguments.
2016-12-07 01:16:36 +00:00
|
|
|
return x + y;
|
2011-09-08 02:10:26 +00:00
|
|
|
# else
|
|
|
|
/* For POWER6 and below isinf/isnan triggers LHS and PLT calls are
|
|
|
|
* costly (especially for POWER6). */
|
|
|
|
# define GET_TW0_HIGH_WORD(d1,d2,i1,i2) \
|
|
|
|
do { \
|
|
|
|
ieee_double_shape_type gh_u1; \
|
|
|
|
ieee_double_shape_type gh_u2; \
|
|
|
|
gh_u1.value = (d1); \
|
|
|
|
gh_u2.value = (d2); \
|
2013-05-17 13:12:16 +00:00
|
|
|
(i1) = gh_u1.parts.msw & 0x7fffffff; \
|
|
|
|
(i2) = gh_u2.parts.msw & 0x7fffffff; \
|
2011-09-08 02:10:26 +00:00
|
|
|
} while (0)
|
|
|
|
|
|
|
|
# define TEST_INF_NAN(x, y) \
|
|
|
|
do { \
|
2013-05-17 13:12:16 +00:00
|
|
|
uint32_t hx, hy; \
|
2011-09-08 02:10:26 +00:00
|
|
|
GET_TW0_HIGH_WORD(x, y, hx, hy); \
|
|
|
|
if (hy > hx) { \
|
|
|
|
uint32_t ht = hx; hx = hy; hy = ht; \
|
|
|
|
} \
|
|
|
|
if (hx >= 0x7ff00000) { \
|
Fix hypot sNaN handling (bug 20940).
TS 18661-1 generally defines libm functions taking sNaN arguments to
return qNaN and raise "invalid", even for the cases where a
corresponding qNaN argument would not result in a qNaN return. This
includes hypot with one argument being an infinity and the other being
an sNaN. This patch duly fixes hypot implementatations in glibc
(generic and powerpc) to ensure qNaN, computed by arithmetic on the
arguments, is returned in that case.
Various implementations do their checks for infinities and NaNs inline
by manipulating the representations of the arguments. For simplicity,
this patch just uses issignaling to check for sNaN arguments. This
could be inlined like the existing code (with due care about reversed
quiet NaN conventions, for implementations where that is relevant),
but given that all these checks are in cases where it's already known
at least one argument is not finite, which should be the uncommon
case, that doesn't seem worthwhile unless performance issues are
observed in practice.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #20940]
* sysdeps/ieee754/dbl-64/e_hypot.c (__ieee754_hypot): Do not
return Inf for arguments Inf and sNaN.
* sysdeps/ieee754/flt-32/e_hypotf.c (__ieee754_hypotf): Likewise.
* sysdeps/ieee754/ldbl-128/e_hypotl.c (__ieee754_hypotl):
Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_hypotl.c (__ieee754_hypotl):
Likewise.
* sysdeps/ieee754/ldbl-96/e_hypotl.c (__ieee754_hypotl): Likewise.
* sysdeps/powerpc/fpu/e_hypot.c (TEST_INF_NAN): Do not return Inf
for arguments Inf and sNaN. When returning a NaN, compute it by
arithmetic on the arguments.
* sysdeps/powerpc/fpu/e_hypotf.c (TEST_INF_NAN): Likewise.
* math/libm-test.inc (pow_test_data): Add tests of sNaN arguments.
2016-12-07 01:16:36 +00:00
|
|
|
if ((hx == 0x7ff00000 || hy == 0x7ff00000) \
|
|
|
|
&& !issignaling (x) && !issignaling (y)) \
|
2011-09-08 02:10:26 +00:00
|
|
|
return INFINITY; \
|
Fix hypot sNaN handling (bug 20940).
TS 18661-1 generally defines libm functions taking sNaN arguments to
return qNaN and raise "invalid", even for the cases where a
corresponding qNaN argument would not result in a qNaN return. This
includes hypot with one argument being an infinity and the other being
an sNaN. This patch duly fixes hypot implementatations in glibc
(generic and powerpc) to ensure qNaN, computed by arithmetic on the
arguments, is returned in that case.
Various implementations do their checks for infinities and NaNs inline
by manipulating the representations of the arguments. For simplicity,
this patch just uses issignaling to check for sNaN arguments. This
could be inlined like the existing code (with due care about reversed
quiet NaN conventions, for implementations where that is relevant),
but given that all these checks are in cases where it's already known
at least one argument is not finite, which should be the uncommon
case, that doesn't seem worthwhile unless performance issues are
observed in practice.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #20940]
* sysdeps/ieee754/dbl-64/e_hypot.c (__ieee754_hypot): Do not
return Inf for arguments Inf and sNaN.
* sysdeps/ieee754/flt-32/e_hypotf.c (__ieee754_hypotf): Likewise.
* sysdeps/ieee754/ldbl-128/e_hypotl.c (__ieee754_hypotl):
Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_hypotl.c (__ieee754_hypotl):
Likewise.
* sysdeps/ieee754/ldbl-96/e_hypotl.c (__ieee754_hypotl): Likewise.
* sysdeps/powerpc/fpu/e_hypot.c (TEST_INF_NAN): Do not return Inf
for arguments Inf and sNaN. When returning a NaN, compute it by
arithmetic on the arguments.
* sysdeps/powerpc/fpu/e_hypotf.c (TEST_INF_NAN): Likewise.
* math/libm-test.inc (pow_test_data): Add tests of sNaN arguments.
2016-12-07 01:16:36 +00:00
|
|
|
return x + y; \
|
2011-09-08 02:10:26 +00:00
|
|
|
} \
|
|
|
|
} while (0)
|
|
|
|
|
|
|
|
#endif
|
|
|
|
|
|
|
|
|
|
|
|
double
|
|
|
|
__ieee754_hypot (double x, double y)
|
|
|
|
{
|
|
|
|
x = fabs (x);
|
|
|
|
y = fabs (y);
|
|
|
|
|
|
|
|
TEST_INF_NAN (x, y);
|
|
|
|
|
|
|
|
if (y > x)
|
|
|
|
{
|
|
|
|
double t = x;
|
|
|
|
x = y;
|
|
|
|
y = t;
|
|
|
|
}
|
2013-05-06 19:40:17 +00:00
|
|
|
if (y == 0.0)
|
|
|
|
return x;
|
|
|
|
/* if y is higher enough, y * 2^60 might overflow. The tests if
|
|
|
|
y >= 1.7976931348623157e+308/2^60 (two60factor) and uses the
|
|
|
|
appropriate check to avoid the overflow exception generation. */
|
|
|
|
if (y > two60factor)
|
2011-09-08 02:10:26 +00:00
|
|
|
{
|
2013-05-06 19:40:17 +00:00
|
|
|
if ((x / y) > two60)
|
|
|
|
return x + y;
|
|
|
|
}
|
|
|
|
else
|
|
|
|
{
|
|
|
|
if (x > (y * two60))
|
|
|
|
return x + y;
|
2011-09-08 02:10:26 +00:00
|
|
|
}
|
|
|
|
if (x > two500)
|
|
|
|
{
|
|
|
|
x *= twoM600;
|
|
|
|
y *= twoM600;
|
2011-09-15 09:16:03 +00:00
|
|
|
return __ieee754_sqrt (x * x + y * y) / twoM600;
|
2011-09-08 02:10:26 +00:00
|
|
|
}
|
|
|
|
if (y < twoM500)
|
|
|
|
{
|
|
|
|
if (y <= pdnum)
|
|
|
|
{
|
|
|
|
x *= two1022;
|
|
|
|
y *= two1022;
|
Fix hypot missing underflows (bug 18803).
Similar to various other bugs in this area, hypot functions can fail
to raise the underflow exception when the result is tiny and inexact
but one or more low bits of the intermediate result that is scaled
down (or, in the i386 case, converted from a wider evaluation format)
are zero. This patch forces the exception in a similar way to
previous fixes.
Note that this issue cannot arise for implementations of hypotf using
double (or wider) for intermediate evaluation (if hypotf should
underflow, that means the double square root is being computed of some
number of the form N*2^-298, for 0 < N < 2^46, which is exactly
represented as a double, and whatever the rounding mode such a square
root cannot have a mantissa with all zeroes after the initial 23
bits). Thus no changes are made to hypotf implementations in this
patch, only to hypot and hypotl.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #18803]
* sysdeps/i386/fpu/e_hypot.S: Use DEFINE_DBL_MIN.
(MO): New macro.
(__ieee754_hypot) [PIC]: Load PIC register.
(__ieee754_hypot): Use DBL_NARROW_EVAL_UFLOW_NONNEG instead of
DBL_NARROW_EVAL.
* sysdeps/ieee754/dbl-64/e_hypot.c (__ieee754_hypot): Use
math_check_force_underflow_nonneg in case where result might be
tiny.
* sysdeps/ieee754/ldbl-128/e_hypotl.c (__ieee754_hypotl):
Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_hypotl.c (__ieee754_hypotl):
Likewise.
* sysdeps/ieee754/ldbl-96/e_hypotl.c (__ieee754_hypotl): Likewise.
* sysdeps/powerpc/fpu/e_hypot.c (__ieee754_hypot): Likewise.
* math/auto-libm-test-in: Add more tests of hypot.
* math/auto-libm-test-out: Regenerated.
2015-09-24 23:43:57 +00:00
|
|
|
double ret = __ieee754_sqrt (x * x + y * y) / two1022;
|
|
|
|
math_check_force_underflow_nonneg (ret);
|
|
|
|
return ret;
|
2011-09-08 02:10:26 +00:00
|
|
|
}
|
|
|
|
else
|
|
|
|
{
|
|
|
|
x *= two600;
|
|
|
|
y *= two600;
|
2011-09-15 09:16:03 +00:00
|
|
|
return __ieee754_sqrt (x * x + y * y) / two600;
|
2011-09-08 02:10:26 +00:00
|
|
|
}
|
|
|
|
}
|
2011-09-15 09:16:03 +00:00
|
|
|
return __ieee754_sqrt (x * x + y * y);
|
2011-09-08 02:10:26 +00:00
|
|
|
}
|
2011-10-12 15:27:51 +00:00
|
|
|
strong_alias (__ieee754_hypot, __hypot_finite)
|