diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/hypot.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/hypot.c | 28 |
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 (;;) { |