1/* mpfr_hypot -- Euclidean distance 2 3Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#define MPFR_NEED_LONGLONG_H 24#include "mpfr-impl.h" 25 26/* The computation of hypot of x and y is done by * 27 * hypot(x,y)= sqrt(x^2+y^2) = z */ 28 29int 30mpfr_hypot (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode) 31{ 32 int inexact, exact; 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; 37 mpfr_exp_t Ex, sh; 38 mpfr_uexp_t diff_exp; 39 40 MPFR_SAVE_EXPO_DECL (expo); 41 MPFR_ZIV_DECL (loop); 42 MPFR_BLOCK_DECL (flags); 43 44 MPFR_LOG_FUNC 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)); 50 51 /* particular cases */ 52 if (MPFR_ARE_SINGULAR (x, y)) 53 { 54 if (MPFR_IS_INF (x) || MPFR_IS_INF (y)) 55 { 56 /* Return +inf, even when the other number is NaN. */ 57 MPFR_SET_INF (z); 58 MPFR_SET_POS (z); 59 MPFR_RET (0); 60 } 61 else if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y)) 62 { 63 MPFR_SET_NAN (z); 64 MPFR_RET_NAN; 65 } 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); 70 } 71 72 if (mpfr_cmpabs (x, y) < 0) 73 { 74 mpfr_srcptr u; 75 u = x; 76 x = y; 77 y = u; 78 } 79 80 /* now |x| >= |y| */ 81 82 Ex = MPFR_GET_EXP (x); 83 diff_exp = (mpfr_uexp_t) Ex - MPFR_GET_EXP (y); 84 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 */ 90 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) */ 95 { 96 if (MPFR_UNLIKELY (rnd_mode == MPFR_RNDU)) 97 { 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) 101 mpfr_nexttoinf (z); 102 MPFR_RET (1); 103 } 104 else /* MPFR_RNDZ, MPFR_RNDD, MPFR_RNDN */ 105 { 106 if (MPFR_LIKELY (Nz >= N)) 107 { 108 mpfr_abs (z, x, rnd_mode); /* exact */ 109 MPFR_RET (-1); 110 } 111 else 112 { 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, 116 goto addoneulp, 117 if (MPFR_UNLIKELY (++ MPFR_EXP (z) > 118 __gmpfr_emax)) 119 return mpfr_overflow (z, rnd_mode, 1); 120 ); 121 122 if (MPFR_UNLIKELY (inexact == 0)) 123 inexact = -1; 124 MPFR_RET (inexact); 125 } 126 } 127 } 128 129 /* General case */ 130 131 N = MAX (MPFR_PREC (x), MPFR_PREC (y)); 132 133 /* working precision */ 134 Nt = Nz + MPFR_INT_CEIL_LOG2 (Nz) + 4; 135 136 mpfr_init2 (t, Nt); 137 mpfr_init2 (te, Nt); 138 mpfr_init2 (ti, Nt); 139 140 MPFR_SAVE_EXPO_MARK (expo); 141 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; 147 148 MPFR_ZIV_INIT (loop, Nt); 149 for (;;) 150 { 151 mpfr_prec_t err; 152 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); 159 160 err = Nt < N ? 4 : 2; 161 if (MPFR_LIKELY (exact == 0 162 || MPFR_CAN_ROUND (t, Nt-err, Nz, rnd_mode))) 163 break; 164 165 MPFR_ZIV_NEXT (loop, Nt); 166 mpfr_set_prec (t, Nt); 167 mpfr_set_prec (te, Nt); 168 mpfr_set_prec (ti, Nt); 169 } 170 MPFR_ZIV_FREE (loop); 171 172 MPFR_BLOCK (flags, inexact = mpfr_div_2si (z, t, sh, rnd_mode)); 173 MPFR_ASSERTD (exact == 0 || inexact != 0); 174 175 mpfr_clear (t); 176 mpfr_clear (ti); 177 mpfr_clear (te); 178 179 /* 180 exact inexact 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 185 */ 186 187 MPFR_SAVE_EXPO_FREE (expo); 188 189 if (MPFR_OVERFLOW (flags)) 190 mpfr_set_overflow (); 191 /* hypot(x,y) >= |x|, thus underflow is not possible. */ 192 193 return mpfr_check_range (z, inexact, rnd_mode); 194} 195