1/* mpfr_hypot -- Euclidean distance 2 3Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao 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 /* particular cases */ 45 if (MPFR_ARE_SINGULAR (x, y)) 46 { 47 if (MPFR_IS_INF (x) || MPFR_IS_INF (y)) 48 { 49 /* Return +inf, even when the other number is NaN. */ 50 MPFR_SET_INF (z); 51 MPFR_SET_POS (z); 52 MPFR_RET (0); 53 } 54 else if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y)) 55 { 56 MPFR_SET_NAN (z); 57 MPFR_RET_NAN; 58 } 59 else if (MPFR_IS_ZERO (x)) 60 return mpfr_abs (z, y, rnd_mode); 61 else /* y is necessarily 0 */ 62 return mpfr_abs (z, x, rnd_mode); 63 } 64 65 if (mpfr_cmpabs (x, y) < 0) 66 { 67 mpfr_srcptr u; 68 u = x; 69 x = y; 70 y = u; 71 } 72 73 /* now |x| >= |y| */ 74 75 Ex = MPFR_GET_EXP (x); 76 diff_exp = (mpfr_uexp_t) Ex - MPFR_GET_EXP (y); 77 78 N = MPFR_PREC (x); /* Precision of input variable */ 79 Nz = MPFR_PREC (z); /* Precision of output variable */ 80 threshold = (MAX (N, Nz) + (rnd_mode == MPFR_RNDN ? 1 : 0)) << 1; 81 if (rnd_mode == MPFR_RNDA) 82 rnd_mode = MPFR_RNDU; /* since the result is positive, RNDA = RNDU */ 83 84 /* Is |x| a suitable approximation to the precision Nz ? 85 (see algorithms.tex for explanations) */ 86 if (diff_exp > threshold) 87 /* result is |x| or |x|+ulp(|x|,Nz) */ 88 { 89 if (MPFR_UNLIKELY (rnd_mode == MPFR_RNDU)) 90 { 91 /* If z > abs(x), then it was already rounded up; otherwise 92 z = abs(x), and we need to add one ulp due to y. */ 93 if (mpfr_abs (z, x, rnd_mode) == 0) 94 mpfr_nexttoinf (z); 95 MPFR_RET (1); 96 } 97 else /* MPFR_RNDZ, MPFR_RNDD, MPFR_RNDN */ 98 { 99 if (MPFR_LIKELY (Nz >= N)) 100 { 101 mpfr_abs (z, x, rnd_mode); /* exact */ 102 MPFR_RET (-1); 103 } 104 else 105 { 106 MPFR_SET_EXP (z, Ex); 107 MPFR_SET_SIGN (z, 1); 108 MPFR_RNDRAW_GEN (inexact, z, MPFR_MANT (x), N, rnd_mode, 1, 109 goto addoneulp, 110 if (MPFR_UNLIKELY (++ MPFR_EXP (z) > 111 __gmpfr_emax)) 112 return mpfr_overflow (z, rnd_mode, 1); 113 ); 114 115 if (MPFR_UNLIKELY (inexact == 0)) 116 inexact = -1; 117 MPFR_RET (inexact); 118 } 119 } 120 } 121 122 /* General case */ 123 124 N = MAX (MPFR_PREC (x), MPFR_PREC (y)); 125 126 /* working precision */ 127 Nt = Nz + MPFR_INT_CEIL_LOG2 (Nz) + 4; 128 129 mpfr_init2 (t, Nt); 130 mpfr_init2 (te, Nt); 131 mpfr_init2 (ti, Nt); 132 133 MPFR_SAVE_EXPO_MARK (expo); 134 135 /* Scale x and y to avoid overflow/underflow in x^2 and overflow in y^2 136 (as |x| >= |y|). The scaling of y can underflow only when the target 137 precision is huge, otherwise the case would already have been handled 138 by the diff_exp > threshold code. */ 139 sh = mpfr_get_emax () / 2 - Ex - 1; 140 141 MPFR_ZIV_INIT (loop, Nt); 142 for (;;) 143 { 144 mpfr_prec_t err; 145 146 exact = mpfr_mul_2si (te, x, sh, MPFR_RNDZ); 147 exact |= mpfr_mul_2si (ti, y, sh, MPFR_RNDZ); 148 exact |= mpfr_sqr (te, te, MPFR_RNDZ); 149 /* Use fma in order to avoid underflow when diff_exp<=MPFR_EMAX_MAX-2 */ 150 exact |= mpfr_fma (t, ti, ti, te, MPFR_RNDZ); 151 exact |= mpfr_sqrt (t, t, MPFR_RNDZ); 152 153 err = Nt < N ? 4 : 2; 154 if (MPFR_LIKELY (exact == 0 155 || MPFR_CAN_ROUND (t, Nt-err, Nz, rnd_mode))) 156 break; 157 158 MPFR_ZIV_NEXT (loop, Nt); 159 mpfr_set_prec (t, Nt); 160 mpfr_set_prec (te, Nt); 161 mpfr_set_prec (ti, Nt); 162 } 163 MPFR_ZIV_FREE (loop); 164 165 MPFR_BLOCK (flags, inexact = mpfr_div_2si (z, t, sh, rnd_mode)); 166 MPFR_ASSERTD (exact == 0 || inexact != 0); 167 168 mpfr_clear (t); 169 mpfr_clear (ti); 170 mpfr_clear (te); 171 172 /* 173 exact inexact 174 0 0 result is exact, ternary flag is 0 175 0 non zero t is exact, ternary flag given by inexact 176 1 0 impossible (see above) 177 1 non zero ternary flag given by inexact 178 */ 179 180 MPFR_SAVE_EXPO_FREE (expo); 181 182 if (MPFR_OVERFLOW (flags)) 183 mpfr_set_overflow (); 184 /* hypot(x,y) >= |x|, thus underflow is not possible. */ 185 186 return mpfr_check_range (z, inexact, rnd_mode); 187} 188