summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/hypot.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/hypot.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/hypot.c28
1 files changed, 14 insertions, 14 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/hypot.c b/Build/source/libs/mpfr/mpfr-src/src/hypot.c
index 44cb6b3e357..d219d6f4515 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/hypot.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/hypot.c
@@ -1,6 +1,6 @@
/* mpfr_hypot -- Euclidean distance
-Copyright 2001-2019 Free Software Foundation, Inc.
+Copyright 2001-2020 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
@@ -29,7 +29,8 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
int
mpfr_hypot (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
{
- int inexact, exact;
+ int inexact;
+ unsigned int exact; /* Warning: 0 will mean "exact" */
mpfr_t t, te, ti; /* auxiliary variables */
mpfr_prec_t N, Nz; /* size variables */
mpfr_prec_t Nt; /* precision of the intermediary variable */
@@ -69,6 +70,8 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
return mpfr_abs (z, x, rnd_mode);
}
+ /* TODO: It may be sufficient to just compare the exponents.
+ The error analysis would need to be updated. */
if (mpfr_cmpabs (x, y) < 0)
{
mpfr_srcptr u;
@@ -155,20 +158,17 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
round toward zero. Using a larger sh wouldn't guarantee an absence
of overflow. Note that the scaling of y can underflow only when the
target precision is huge, otherwise the case would already have been
- handled by the diff_exp > threshold code.
- FIXME: Friedland in "Algorithm 312: Absolute Value and Square Root of a
- Complex Number" (Communications of the ACM, 1967) avoids overflow by
- computing |x|*sqrt(1+(y/x)^2) if |x| >= |y|, and |y|*sqrt(1+(x/y)^2)
- otherwise.
- [VL] This trick (which is a scaling by a non-power of 2, thus doesn't
- really bring new behavior w.r.t. overflow/underflow exceptions) may be
- useful for hardware floating-point formats because a whole power-of-2
- scaling code is likely to take more time than the additional division,
- but in the context of multiple-precision, I doubt that it is a good
- idea. Ideally scaling by a power of 2 could be done in a constant time,
- e.g. with MPFR_ALIAS; but one needs to be very careful... */
+ handled by the diff_exp > threshold code; but this case is avoided
+ thanks to a FMA (this problem is transferred to the FMA code). */
sh = (mpfr_get_emax () - 1) / 2 - Ex;
+ /* FIXME: ti is subject to underflow. Solution: x and y could be
+ aliased with MPFR_ALIAS, and if need be, the aliases be pre-scaled
+ exactly as UBF, so that x^2 + y^2 is in range. Then call mpfr_fmma
+ and the square root, and scale the result. The error analysis would
+ be simpler.
+ Note: mpfr_fmma is currently not optimized. */
+
MPFR_ZIV_INIT (loop, Nt);
for (;;)
{