1 /* mpfr_hypot -- Euclidean distance
3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 /* The computation of hypot of x and y is done by *
27 * hypot(x,y)= sqrt(x^2+y^2) = z */
30 mpfr_hypot (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
33 mpfr_t t, te, ti; /* auxiliary variables */
34 mpfr_prec_t N, Nz; /* size variables */
35 mpfr_prec_t Nt; /* precision of the intermediary variable */
36 mpfr_prec_t threshold;
40 MPFR_SAVE_EXPO_DECL (expo);
42 MPFR_BLOCK_DECL (flags);
45 (("x[%Pu]=%.*Rg y[%Pu]=%.*Rg rnd=%d",
46 mpfr_get_prec (x), mpfr_log_prec, x,
47 mpfr_get_prec (y), mpfr_log_prec, y, rnd_mode),
48 ("z[%Pu]=%.*Rg inexact=%d",
49 mpfr_get_prec (z), mpfr_log_prec, z, inexact));
51 /* particular cases */
52 if (MPFR_ARE_SINGULAR (x, y))
54 if (MPFR_IS_INF (x) || MPFR_IS_INF (y))
56 /* Return +inf, even when the other number is NaN. */
61 else if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y))
66 else if (MPFR_IS_ZERO (x))
67 return mpfr_abs (z, y, rnd_mode);
68 else /* y is necessarily 0 */
69 return mpfr_abs (z, x, rnd_mode);
72 if (mpfr_cmpabs (x, y) < 0)
82 Ex = MPFR_GET_EXP (x);
83 diff_exp = (mpfr_uexp_t) Ex - MPFR_GET_EXP (y);
85 N = MPFR_PREC (x); /* Precision of input variable */
86 Nz = MPFR_PREC (z); /* Precision of output variable */
87 threshold = (MAX (N, Nz) + (rnd_mode == MPFR_RNDN ? 1 : 0)) << 1;
88 if (rnd_mode == MPFR_RNDA)
89 rnd_mode = MPFR_RNDU; /* since the result is positive, RNDA = RNDU */
91 /* Is |x| a suitable approximation to the precision Nz ?
92 (see algorithms.tex for explanations) */
93 if (diff_exp > threshold)
94 /* result is |x| or |x|+ulp(|x|,Nz) */
96 if (MPFR_UNLIKELY (rnd_mode == MPFR_RNDU))
98 /* If z > abs(x), then it was already rounded up; otherwise
99 z = abs(x), and we need to add one ulp due to y. */
100 if (mpfr_abs (z, x, rnd_mode) == 0)
104 else /* MPFR_RNDZ, MPFR_RNDD, MPFR_RNDN */
106 if (MPFR_LIKELY (Nz >= N))
108 mpfr_abs (z, x, rnd_mode); /* exact */
113 MPFR_SET_EXP (z, Ex);
114 MPFR_SET_SIGN (z, 1);
115 MPFR_RNDRAW_GEN (inexact, z, MPFR_MANT (x), N, rnd_mode, 1,
117 if (MPFR_UNLIKELY (++ MPFR_EXP (z) >
119 return mpfr_overflow (z, rnd_mode, 1);
122 if (MPFR_UNLIKELY (inexact == 0))
131 N = MAX (MPFR_PREC (x), MPFR_PREC (y));
133 /* working precision */
134 Nt = Nz + MPFR_INT_CEIL_LOG2 (Nz) + 4;
140 MPFR_SAVE_EXPO_MARK (expo);
142 /* Scale x and y to avoid overflow/underflow in x^2 and overflow in y^2
143 (as |x| >= |y|). The scaling of y can underflow only when the target
144 precision is huge, otherwise the case would already have been handled
145 by the diff_exp > threshold code. */
146 sh = mpfr_get_emax () / 2 - Ex - 1;
148 MPFR_ZIV_INIT (loop, Nt);
153 exact = mpfr_mul_2si (te, x, sh, MPFR_RNDZ);
154 exact |= mpfr_mul_2si (ti, y, sh, MPFR_RNDZ);
155 exact |= mpfr_sqr (te, te, MPFR_RNDZ);
156 /* Use fma in order to avoid underflow when diff_exp<=MPFR_EMAX_MAX-2 */
157 exact |= mpfr_fma (t, ti, ti, te, MPFR_RNDZ);
158 exact |= mpfr_sqrt (t, t, MPFR_RNDZ);
160 err = Nt < N ? 4 : 2;
161 if (MPFR_LIKELY (exact == 0
162 || MPFR_CAN_ROUND (t, Nt-err, Nz, rnd_mode)))
165 MPFR_ZIV_NEXT (loop, Nt);
166 mpfr_set_prec (t, Nt);
167 mpfr_set_prec (te, Nt);
168 mpfr_set_prec (ti, Nt);
170 MPFR_ZIV_FREE (loop);
172 MPFR_BLOCK (flags, inexact = mpfr_div_2si (z, t, sh, rnd_mode));
173 MPFR_ASSERTD (exact == 0 || inexact != 0);
181 0 0 result is exact, ternary flag is 0
182 0 non zero t is exact, ternary flag given by inexact
183 1 0 impossible (see above)
184 1 non zero ternary flag given by inexact
187 MPFR_SAVE_EXPO_FREE (expo);
189 if (MPFR_OVERFLOW (flags))
190 mpfr_set_overflow ();
191 /* hypot(x,y) >= |x|, thus underflow is not possible. */
193 return mpfr_check_range (z, inexact, rnd_mode);