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
);