Merged changesets concerning mpfr_hypot, from trunk: 5291, 5292, 5297,

5299, 5301, 5302, 5308, 5310, 5311, 5319, 5336, 5464, 5465, 5473, 5563,
5564, 5566. Details:
* hypot.c: rewrite, fixing the following bugs:
  - the inexact flag was not always set as required;
  - mpfr_hypot could loop on tiny values due to internal underflow;
  - for generic overflows, the overflow flag was lost (note: this
    bug is independent of the bug fixed in r5546 since the overflow
    does not come after a can_round).
* tests/thypot.c:
  - changed custom random tests to tgeneric ones;
  - added testcases for the above bugs.
* algorithms.tex: mpfr_hypot proof rewritten.
* BUGS: added a bug concerning mpfr_hypot.

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/2.3@5568 280ebfd0-de03-0410-8827-d642c229c3f4

4 files changed, 310 insertions, 162 deletions

diff --git a/BUGS b/BUGS
index a7bc99c7f..a582c104e 100644
--- a/BUGS
+++ b/BUGS
@@ -37,6 +37,11 @@ Known bugs:
 * Some functions do not use MPFR_SAVE_EXPO_* macros, thus do not behave
   correctly in a reduced exponent range.
 
+* Function hypot gives incorrect result when on the one hand the difference
+  between parameters' exponents is near 2*MPFR_EMAX_MAX and on the other hand
+  the output precision or the precision of the parameter with greatest
+  absolute value is greater than 2*MPFR_EMAX_MAX-4.
+
 Potential bugs:
 
 * Possible incorrect results due to internal underflow, which can lead to
diff --git a/algorithms.tex b/algorithms.tex
index 9567fb01f..125078a6c 100644
--- a/algorithms.tex
+++ b/algorithms.tex
@@ -2168,20 +2168,121 @@ The \texttt{mpfr\_hypot} function implements the euclidean distance function:
 \[
 \textnormal{hypot} (x,y) = \sqrt{x^2+y^2}.
 \]
-The algorithm used is as follows, where all roundings are done towards zero:
-\begin{eqnarray}\nonumber
-u&\leftarrow&\circ(x^2)\\\nonumber
-v&\leftarrow&\circ(y^2)\\\nonumber
-w&\leftarrow&\circ(u+v)\\\nonumber
-s&\leftarrow&\circ(\sqrt{w})
-\end{eqnarray}
-Since the inputs $x$ and $y$ are exact, we have $|u-x^2| \leq \ulp(u)$
-and $|v-y^2| \leq \ulp(v)$;
-then $|w-(x^2+y^2)| \leq \ulp(w) + |u-x^2| + |v-y^2| \leq 3 \ulp(w)$
-since $w$ is greater or equal to $u$ and $v$.
-For the last step $s \leftarrow \circ(\sqrt{w})$, we use the last formula
-from~\textsection\ref{generic:sqrt}, which gives
-$|s-\sqrt{x^2+y^2}| \leq 4 \ulp(s)$.
+If one of the variable is zero, then hypot is computed using the absolute
+value of the other variable.  Assume that $0 < y \leq x$. Using the first
+degree Taylor polynomial, we have
+\[
+0 < \sqrt{x^2+y^2}-x < \frac{y^2}{2x}.
+\]
+
+Let $p_x, p_y$ the precisions of input variables $x$ and $y$ respectively,
+$p_z$ the output precision and $z=\circ_{p_z}(\sqrt{x^2+y^2})$ the expected
+result.  Let us assume, as it is the case in MPFR, that the minimal and
+maximal acceptable exponents (respectively $e_{min}$ and $e_{max}$) verify $2
+< e_{max}$ and $e_{max} = -e_{min}$.
+
+When rounding to nearest, if $p_x \leq p_z$ and $\frac{p_z+1}{2} < \Exp(x) -
+\Exp(y)$, we have $\frac{y^2}{2x} < \frac{1}{2}\ulp_{p_z}(x)$ ; if $p_z <
+p_x$, the condition $\frac{p_x+1}{2} < \Exp(x) - \Exp(y)$ ensures that
+$\frac{y^2}{2x} < \frac{1}{2} \ulp_{p_x}(x)$.  In both cases, these
+inequalities show that $z=\N_{p_z}(x)$, except that tie case is rounded
+towards plus infinity since hypot($x$,$y$) is greater than but not equal to
+$x$.
+
+With other rounding modes, the conditions $p_z/2 < \Exp(x) - \Exp(y)$ if $p_x
+\leq p_z$, and $p_x/2 < \Exp(x) - \Exp(y)$ if $p_z < p_x$ mean in a similar
+way that $z=\circ_{p_z}(x)$, except that we need to add one ulp to the result
+when rounding towards plus infinity and $x$ is exactly representable with
+$p_z$ bits of precision.
+
+When none of the above conditions is satisfied and when $\Exp(x) - \Exp(y)
+\leq e_{max} - 2$, we use the following algorithm:
+
+\begin{quote}
+  Algorithm {\tt hypot\_1}\\
+  Input: $x$ and $y$ with $|y| \leq |x|$,
+  $p$ the working precision with $p \geq p_z$.\\
+  Output: $\sqrt{x^2+y^2}$ with $\left\{
+  \begin{array}{l} 
+    p-4 \textnormal{ bits of precision if } p <\max(p_x, p_y),\\
+    p-2 \textnormal{ bits of precision if } \max(p_x, p_y) \leq p.
+  \end{array}\right.$\\
+  $s \leftarrow \lfloor e_{max}/2\rfloor - \Exp(x) -1$\\
+  $x_s \leftarrow \Z(x\times 2^s)$\\
+  $y_s \leftarrow \Z(y\times 2^s)$\\
+  $u \leftarrow \Z(x_s^2)$\\
+  $v \leftarrow \Z(y_s^2)$\\
+  $w \leftarrow \Z(u+v)$\\
+  $t \leftarrow \Z(\sqrt{w})$\\
+  $z \leftarrow \Z(t/2^s)$
+\end{quote}
+
+In order to avoid undue overflow during computation, we shift inputs'
+exponents by $s = \lfloor\frac{e_{max}}{2}\rfloor -1 -\Exp(x)$ before
+computing squares and shift back the output's exponent by $-s$ using the fact
+that $\sqrt{(x.2^s)^2+(y.2^s)^2}/2^s = \sqrt{x^2+y^2}$. We show below that no
+overflow nor underflow goes on.
+
+We check first that the exponent shift do not cause overflow and, in the same
+time, that the squares of the shifted inputs do never overflow. For $x$, we
+have $\Exp(x) + s = \lfloor e_{max}/2\rfloor - 1$, so $\Exp(x_s^2) < e_{max} -
+1$ and neither $x_s$ nor $x_s^2$ overflows. For $y$, note that we have $y_s
+\leq x_s$ because $y \leq x$, thus $y_s$ and $y_s^2$ do not overflow.
+
+Secondly, let us see that the exponent shift do not cause underflow. For $x$,
+we know that $0 \leq \Exp(x) + s$, thus neither $x_s$ nor $x_s^2$
+underflows. For $y$, the condition $\Exp(x) - \Exp(y) \leq e_{max} - 2$
+ensures that $-e_{max}/2 \leq \Exp(y) + s$ which shows that $y_s$ and its
+square do not underflow.
+
+Thirdly, the addition does not overflow because $u + v < 2 x_s^2$ and it was
+shown above that $x_s^2 < 2^{e_{max} - 1}$. It cannot underflow because both
+operands are positive.
+
+Fourthly, as $x_s < t$, the square root does not underflow. Due to the
+exponent shift, we have $1 \leq x_s$, then $w$ is greater than 1 and thus
+greater than its square root $t$, so the square root does overflow.
+
+Finally, let us show that the back shift do not raise underflow nor overflow
+unless the exact result is greater than or equal to $2^{e_{max}}$. Because no
+underflow has occured so far $\Exp(x) \leq \Exp(t) - s$ which shows that it
+does not underflow. And all roundings being towards zero we have $z \leq
+\sqrt{x^2 + y^2}$, so if $2^{e_{max}} \leq z$, then the exact value is also
+greater than or equal to $2^{e_{max}}$.
+
+Let us analyse now the error of the algorithm hypot\_1:\\
+\begin{tabular}{l p{2em} l c r | r}
+  & & & & $p < \min(p_x, p_y)$ & $\max(p_x, p_y) \leq p$\\
+  $x_s \leftarrow \Z(x\times 2^s)$ & & $\error(x_s)$ & $\leq$ &
+  1 $\ulp(x_s)$ &
+  exact\\
+  $y_s \leftarrow \Z(y\times 2^s)$ & & $\error(y_s)$ & $\leq$ &
+  1 $\ulp(y_s)$ &
+  exact\\
+  $u \leftarrow \Z(x_s^2)$ & & $\error(u)$ & $\leq$ &
+  6 $\ulp(u)$ &
+  1 $\ulp(u)$\\
+  $v \leftarrow \Z(y_s^2)$ & & $\error(v)$ & $\leq$ &
+  6 $\ulp(v)$ &
+  1 $\ulp(v)$\\
+  $w \leftarrow \Z(u+v)$ & & $\error(w)$ & $\leq$ &
+  13 $\ulp(w)$ &
+  3 $\ulp(w)$\\
+  $t \leftarrow \Z(\sqrt{w})$ & & $\error(t)$ & $\leq$ &
+  14 $\ulp(t)$ &
+  4 $\ulp(t)$\\
+  $z \leftarrow \Z(t/2^s)$& & exact.\\
+\end{tabular}\\
+And in the intermediate case, if $\min(p_x,p_y) \leq p < \max(p_x,p_y)$, we
+have\\
+\begin{tabular}{l p{2em} l c r l}
+  $w \leftarrow \Z(u+v)$ & &
+  $\error(w)$ & $\leq$ & 8 $\ulp(w)$\\
+  $t \leftarrow \Z(\sqrt{w})$ & &
+  $\error(t)$ & $\leq$ & 9 $\ulp(t)$.\\
+\end{tabular}\\
+Thus, 2 bits of precision are lost when $\max(p_x, p_y) \leq p$ and 4 bits
+when $p$ does not verify this relation.
 
 \subsection{The floating multiply-add function}
 
diff --git a/hypot.c b/hypot.c
index f9a00f8d7..782be4e8c 100644
--- a/hypot.c
+++ b/hypot.c
@@ -31,12 +31,15 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
 {
   int inexact, exact;
   mpfr_t t, te, ti; /* auxiliary variables */
-  mp_prec_t Nx, Ny, Nz; /* size variables */
+  mp_prec_t N, Nz; /* size variables */
   mp_prec_t Nt;   /* precision of the intermediary variable */
-  mp_exp_t Ex, Ey, sh;
+  mp_prec_t threshold;
+  mp_exp_t Ex, sh;
   mp_exp_unsigned_t diff_exp;
+
   MPFR_SAVE_EXPO_DECL (expo);
   MPFR_ZIV_DECL (loop);
+  MPFR_BLOCK_DECL (flags);
 
   /* particular cases */
   if (MPFR_ARE_SINGULAR (x, y))
@@ -71,25 +74,16 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
   /* now |x| >= |y| */
 
   Ex = MPFR_GET_EXP (x);
-  Ey = MPFR_GET_EXP (y);
-  diff_exp = (mp_exp_unsigned_t) Ex - Ey;
+  diff_exp = (mp_exp_unsigned_t) Ex - MPFR_GET_EXP (y);
 
-  Nx = MPFR_PREC (x);   /* Precision of input variable */
+  N = MPFR_PREC (x);   /* Precision of input variable */
   Nz = MPFR_PREC (z);   /* Precision of output variable */
+  threshold =
+    rnd_mode == GMP_RNDN ? (MAX (N, Nz) + 1) <<1 : MAX (N, Nz) <<1;
 
-  /* we have x < 2^Ex thus x^2 < 2^(2*Ex),
-     and ulp(x) = 2^(Ex-Nx) thus ulp(x^2) >= 2^(2*Ex-2*Nx).
-     y does not overlap with the result when
-     x^2+y^2 < (|x| + 1/2*ulp(x,Nz))^2 = x^2 + |x|*ulp(x,Nz) + 1/4*ulp(x,Nz)^2,
-     i.e. a sufficient condition is y^2 < |x|*ulp(x,Nz),
-     or 2^(2*Ey) <= 2^(2*Ex-1-Nz), i.e. 2*diff_exp > Nz.
-     Warning: this is true only for Nx <= Nz, otherwise the trailing bits
-     of x may be already very close to 1/2*ulp(x,Nz)!
-     If Nx > Nz, then we can notice that it is possible to round on Nx bits
-     if 2*diff_exp > Nx (see above as if Nz = Nx), therefore on Nz bits.
-     Hence the condition: 2*diff_exp > MAX(Nz,Nx).
-  */
-  if (diff_exp > MAX (Nz, Nx) / 2)
+  /* Is |x| a suitable approximation to the precision Nz ?
+     (see algorithms.tex for explanations) */
+  if (diff_exp > threshold)
     /* result is |x| or |x|+ulp(|x|,Nz) */
     {
       if (MPFR_UNLIKELY (rnd_mode == GMP_RNDU))
@@ -98,73 +92,76 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
              z = abs(x), and we need to add one ulp due to y. */
           if (mpfr_abs (z, x, rnd_mode) == 0)
             mpfr_nexttoinf (z);
-          return 1;
+          MPFR_RET (1);
         }
       else /* GMP_RNDZ, GMP_RNDD, GMP_RNDN */
         {
-          if (MPFR_LIKELY (Nz >= Nx))
+          if (MPFR_LIKELY (Nz >= N))
             {
               mpfr_abs (z, x, rnd_mode);  /* exact */
-              return -1;
+              MPFR_RET (-1);
             }
           else
             {
               MPFR_SET_EXP (z, Ex);
               MPFR_SET_SIGN (z, 1);
-              MPFR_RNDRAW_GEN (inexact, z, MPFR_MANT (x), Nx, rnd_mode, 1,
+              MPFR_RNDRAW_GEN (inexact, z, MPFR_MANT (x), N, rnd_mode, 1,
                                goto addoneulp,
-                  if (MPFR_UNLIKELY (++MPFR_EXP (z) > __gmpfr_emax))
-                    return mpfr_overflow (z, rnd_mode, 1);
-                              );
-              return inexact ? inexact : -1;
+                               if (MPFR_UNLIKELY (++MPFR_EXP (z) > __gmpfr_emax))
+                                 return mpfr_overflow (z, rnd_mode, 1);
+                               );
+
+              if (MPFR_UNLIKELY (inexact == 0))
+                inexact = -1;
+              MPFR_RET (inexact);
             }
         }
     }
 
   /* General case */
 
-  Ny = MPFR_PREC(y);   /* Precision of input variable */
+  /* FIXME: the following algorithm is correct only for
+     diff_exp <= MPFR_EMAX_MAX - 2 */
+
+  N = MAX (MPFR_PREC (x), MPFR_PREC (y));
 
-  /* compute the working precision -- see algorithms.ps */
-  Nt = MAX (MAX (MAX (Nx, Ny), Nz), 8);
-  /* FIXME: if Nx or Ny are very large with respect to the target precision
-     Nz, this may be overkill! */
-  Nt = Nt + MPFR_INT_CEIL_LOG2 (Nt) + 2;
+  /* working precision */
+  Nt = Nz + MPFR_INT_CEIL_LOG2 (Nz) + 4;
 
-  /* initialise the intermediary variables */
   mpfr_init2 (t, Nt);
   mpfr_init2 (te, Nt);
   mpfr_init2 (ti, Nt);
 
   MPFR_SAVE_EXPO_MARK (expo);
 
-  sh = MAX (0, MIN (Ex, Ey));
+  /* Scale x and y to avoid overflow/underflow in x^2 and y^2. */
+  sh = mpfr_get_emax()/2 -Ex -1;
 
   MPFR_ZIV_INIT (loop, Nt);
   for (;;)
     {
-      /* computations of hypot */
-      mpfr_div_2ui (te, x, sh, GMP_RNDZ); /* exact since Nt >= Nx */
-      mpfr_div_2ui (ti, y, sh, GMP_RNDZ); /* exact since Nt >= Ny */
-      exact = mpfr_mul (te, te, te, GMP_RNDZ);    /* x^2 */
-      exact |= mpfr_mul (ti, ti, ti, GMP_RNDZ);   /* y^2 */
-      exact |= mpfr_add (t, te, ti, GMP_RNDZ);    /* x^2+y^2 */
-      exact |= mpfr_sqrt (t, t, GMP_RNDZ);        /* sqrt(x^2+y^2)*/
+      mp_prec_t err;
 
+      exact = mpfr_mul_2si (te, x, sh, GMP_RNDZ);
+      exact |= mpfr_mul_2si (ti, y, sh, GMP_RNDZ);
+      exact |= mpfr_sqr (te, te, GMP_RNDZ);
+      exact |= mpfr_sqr (ti, ti, GMP_RNDZ);
+      exact |= mpfr_add (t, te, ti, GMP_RNDZ);
+      exact |= mpfr_sqrt (t, t, GMP_RNDZ);
+
+      err = Nt < N ? 4 : 2;
       if (MPFR_LIKELY (exact == 0
-                       || MPFR_CAN_ROUND (t, Nt-2, Nz, rnd_mode)))
+                       || MPFR_CAN_ROUND (t, Nt-err, Nz, rnd_mode)))
         break;
 
-      /* reactualization of the precision */
       MPFR_ZIV_NEXT (loop, Nt);
       mpfr_set_prec (t, Nt);
       mpfr_set_prec (te, Nt);
       mpfr_set_prec (ti, Nt);
-
     }
   MPFR_ZIV_FREE (loop);
-  inexact = mpfr_mul_2ui (z, t, sh, rnd_mode);
 
+  MPFR_BLOCK (flags, inexact = mpfr_div_2si (z, t, sh, rnd_mode));
   MPFR_ASSERTD (exact == 0 || inexact != 0);
 
   mpfr_clear (t);
@@ -172,14 +169,18 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
   mpfr_clear (te);
 
   /*
-       exact  inexact
-        0         0         result is exact, ternary flag is 0
-        0       non zero    t is exact, ternary flag given by inexact
-        1         0         impossible (see above)
-        1       non zero    ternary flag given by inexact
-   */
+    exact  inexact
+    0         0         result is exact, ternary flag is 0
+    0       non zero    t is exact, ternary flag given by inexact
+    1         0         impossible (see above)
+    1       non zero    ternary flag given by inexact
+  */
 
   MPFR_SAVE_EXPO_FREE (expo);
 
+  if (MPFR_OVERFLOW (flags))
+    mpfr_set_overflow ();
+  /* hypot(x,y) >= |x|, thus underflow is not possible. */
+
   return mpfr_check_range (z, inexact, rnd_mode);
 }
diff --git a/tests/thypot.c b/tests/thypot.c
index a13ef030b..69fb90942 100644
--- a/tests/thypot.c
+++ b/tests/thypot.c
@@ -26,7 +26,8 @@ MA 02110-1301, USA. */
 
 #include "mpfr-test.h"
 
-#define TEST_FUNCTION mpfr_hypot
+/* Non-zero when extended exponent range */
+static int ext = 0;
 
 static void
 special (void)
@@ -94,19 +95,34 @@ test_large (void)
   mpfr_set_prec (t, 53);
   mpfr_set_prec (y, 53);
   mpfr_set_str_binary (x, "0.11101100011110000011101000010101010011001101000001100E-1021");
-  mpfr_set_str_binary (t, "0.11111001010011000001110110001101011100001000010010100E-1021");
-  mpfr_hypot (y, x, t, GMP_RNDN);
+  mpfr_set_str_binary (y, "0.11111001010011000001110110001101011100001000010010100E-1021");
+  mpfr_hypot (t, x, y, GMP_RNDN);
+  mpfr_set_str_binary (z, "0.101010111100110111101110111110100110010011001010111E-1020");
+  if (mpfr_cmp (z, t))
+    {
+      printf ("Error in test_large: got\n");
+      mpfr_out_str (stdout, 2, 0, t, GMP_RNDN);
+      printf ("\ninstead of\n");
+      mpfr_out_str (stdout, 2, 0, z, GMP_RNDN);
+      printf ("\n");
+      exit (1);
+    }
 
   mpfr_set_prec (x, 240);
   mpfr_set_prec (y, 22);
+  mpfr_set_prec (z, 2);
   mpfr_set_prec (t, 2);
   mpfr_set_str_binary (x, "0.100111011010010010110100000100000001100010011100110101101111111101011110111011011101010110100101111000111100010100110000100101011110111011100110100110100101110101101100011000001100000001111101110100100100011011011010110111100110010101000111e-7");
-  mpfr_set_str_binary (y, "0.1111000010000011000111e-10");
+  mpfr_set_str_binary (y, "0.1111000010000011000111E-10");
   mpfr_hypot (t, x, y, GMP_RNDN);
-  mpfr_set_str_binary (y, "0.11E-7");
-  if (mpfr_cmp (t, y))
+  mpfr_set_str_binary (z, "0.11E-7");
+  if (mpfr_cmp (z, t))
     {
-      printf ("Error in mpfr_hypot (1)\n");
+      printf ("Error in test_large: got\n");
+      mpfr_out_str (stdout, 2, 0, t, GMP_RNDN);
+      printf ("\ninstead of\n");
+      mpfr_out_str (stdout, 2, 0, z, GMP_RNDN);
+      printf ("\n");
       exit (1);
     }
 
@@ -117,6 +133,50 @@ test_large (void)
 }
 
 static void
+test_small (void)
+{
+  mpfr_t x, y, z1, z2;
+  int inex1, inex2;
+  unsigned int flags;
+
+  /* Test hypot(x,x) with x = 2^(emin-1). Result is x * sqrt(2). */
+  mpfr_inits2 (8, x, y, z1, z2, (mpfr_ptr) 0);
+  mpfr_set_si_2exp (x, 1, mpfr_get_emin () - 1, GMP_RNDN);
+  mpfr_set_si_2exp (y, 1, mpfr_get_emin () - 1, GMP_RNDN);
+  mpfr_set_ui (z1, 2, GMP_RNDN);
+  inex1 = mpfr_sqrt (z1, z1, GMP_RNDN);
+  inex2 = mpfr_mul (z1, z1, x, GMP_RNDN);
+  MPFR_ASSERTN (inex2 == 0);
+  mpfr_clear_flags ();
+  inex2 = mpfr_hypot (z2, x, y, GMP_RNDN);
+  flags = __gmpfr_flags;
+  if (mpfr_cmp (z1, z2) != 0)
+    {
+      printf ("Error in test_small%s\nExpected ",
+              ext ? ", extended exponent range" : "");
+      mpfr_out_str (stdout, 2, 0, z1, GMP_RNDN);
+      printf ("\nGot      ");
+      mpfr_out_str (stdout, 2, 0, z2, GMP_RNDN);
+      printf ("\n");
+      exit (1);
+    }
+  if (! SAME_SIGN (inex1, inex2))
+    {
+      printf ("Bad ternary value in test_small%s\nExpected %d, got %d\n",
+              ext ? ", extended exponent range" : "", inex1, inex2);
+      exit (1);
+    }
+  if (flags != MPFR_FLAGS_INEXACT)
+    {
+      printf ("Bad flags in test_small%s\nExpected %u, got %u\n",
+              ext ? ", extended exponent range" : "",
+              (unsigned int) MPFR_FLAGS_INEXACT, flags);
+      exit (1);
+    }
+  mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0);
+}
+
+static void
 test_large_small (void)
 {
   mpfr_t x, y, z;
@@ -131,7 +191,8 @@ test_large_small (void)
   inexact = mpfr_hypot (z, x, y, GMP_RNDN);
   if (inexact >= 0 || mpfr_cmp (x, z))
     {
-      printf ("Error 1 in test_large_small\n");
+      printf ("Error 1 in test_large_small%s\n",
+              ext ? ", extended exponent range" : "");
       exit (1);
     }
 
@@ -139,7 +200,8 @@ test_large_small (void)
   inexact = mpfr_hypot (z, x, y, GMP_RNDN);
   if (mpfr_cmp (x, z) >= 0)
     {
-      printf ("Error 2 in test_large_small\n");
+      printf ("Error 2 in test_large_small%s\n",
+              ext ? ", extended exponent range" : "");
       printf ("x = ");
       mpfr_out_str (stdout, 2, 0, x, GMP_RNDN);
       printf ("\n");
@@ -159,105 +221,84 @@ test_large_small (void)
   mpfr_clear (z);
 }
 
-int
-main (int argc, char *argv[])
+static void
+check_overflow (void)
 {
-  unsigned int prec, err, yprec, n, p0 = 2, p1 = 100, N = 100;
-  mp_rnd_t rnd;
-  mpfr_t x1, x2, y, z, t;
-  int inexact, compare, compare2;
+  mpfr_t x, y;
+  int inex, r;
 
-  tests_start_mpfr ();
+  mpfr_inits2 (8, x, y, (mpfr_ptr) 0);
+  mpfr_set_ui (x, 1, GMP_RNDN);
+  mpfr_setmax (x, mpfr_get_emax ());
 
-  special ();
-
-  mpfr_init (x1);
-  mpfr_init (x2);
-  mpfr_init (y);
-  mpfr_init (z);
-  mpfr_init (t);
-
-  /* thypot prec - perform one random computation with precision prec */
-  if (argc >= 2)
+  RND_LOOP(r)
     {
-      p0 = p1 = atoi (argv[1]);
-      N = 1;
+      mpfr_clear_overflow ();
+      inex = mpfr_hypot (y, x, x, r);
+      if (!mpfr_overflow_p ())
+        {
+          printf ("No overflow in check_overflow for %s%s\n",
+                  mpfr_print_rnd_mode ((mp_rnd_t) r),
+                  ext ? ", extended exponent range" : "");
+          exit (1);
+        }
+      MPFR_ASSERTN (MPFR_IS_POS (y));
+      if (r == GMP_RNDZ || r == GMP_RNDD)
+        {
+          MPFR_ASSERTN (inex < 0);
+          MPFR_ASSERTN (!mpfr_inf_p (y));
+          mpfr_nexttoinf (y);
+        }
+      else
+        {
+          MPFR_ASSERTN (inex > 0);
+        }
+      MPFR_ASSERTN (mpfr_inf_p (y));
     }
 
-  for (prec = p0; prec <= p1; prec++)
-    {
-      mpfr_set_prec (x1, prec);
-      mpfr_set_prec (x2, prec);
-      mpfr_set_prec (z, prec);
-      mpfr_set_prec (t, prec);
-      yprec = prec + 10;
+  mpfr_clears (x, y, (mpfr_ptr) 0);
+}
 
-      for (n=0; n<N; n++)
-        {
-          mpfr_random(x1);
-          mpfr_random(x2);
-          if (randlimb () % 2)
-            mpfr_neg (x1, x1, GMP_RNDN);
-          if (randlimb () % 2)
-            mpfr_neg (x2, x2, GMP_RNDN);
-          rnd = (mp_rnd_t) RND_RAND ();
-          mpfr_set_prec (y, yprec);
-
-          compare =TEST_FUNCTION (y, x1,x2, rnd);
-          err = (rnd == GMP_RNDN) ? yprec + 1 : yprec;
-          if (mpfr_can_round (y, err, rnd, rnd, prec))
-            {
-              mpfr_set (t, y, rnd);
-              inexact = TEST_FUNCTION (z, x1,x2, rnd);
-              if (mpfr_cmp (t, z))
-                {
-                  printf ("results differ for x1=");
-                  mpfr_out_str (stdout, 2, prec, x1, GMP_RNDN);
-                  printf ("\n and x2=");
-                  mpfr_out_str (stdout, 2, prec, x2, GMP_RNDN);
-                  printf (" \n prec=%u rnd_mode=%s\n", prec,
-                          mpfr_print_rnd_mode (rnd));
-                  printf ("   got ");
-                  mpfr_out_str (stdout, 2, prec, z, GMP_RNDN);
-                  puts ("");
-                  printf ("   expected ");
-                  mpfr_out_str (stdout, 2, prec, t, GMP_RNDN);
-                  puts ("");
-                  printf ("   approximation was ");
-                  mpfr_print_binary (y);
-                  puts ("");
-                  exit (1);
-                }
-              compare2 = mpfr_cmp (t, y);
-              /* if rounding to nearest, cannot know the sign of t - f(x)
-                 because of composed rounding: y = o(f(x)) and t = o(y) */
-              if ((rnd != GMP_RNDN) && (compare * compare2 >= 0))
-                compare = compare + compare2;
-              else
-                compare = inexact; /* cannot determine sign(t-f(x)) */
-              if (((inexact == 0) && (compare != 0)) ||
-                  ((inexact > 0) && (compare <= 0)) ||
-                  ((inexact < 0) && (compare >= 0)))
-                {
-                  printf ("Wrong inexact flag for rnd=%s: expected %d, got %d"
-                          "\n", mpfr_print_rnd_mode (rnd), compare, inexact);
-                  printf ("x1="); mpfr_print_binary (x1); puts ("");
-                  printf ("x2="); mpfr_print_binary (x2); puts ("");
-                  printf ("t="); mpfr_print_binary (t); puts ("");
-                  exit (1);
-                }
-            }
-        }
+#define TWO_ARGS
+#define TEST_FUNCTION mpfr_hypot
+#include "tgeneric.c"
+
+static void
+alltst (void)
+{
+  mp_exp_t emin, emax;
+
+  ext = 0;
+  test_small ();
+  test_large_small ();
+  check_overflow ();
+
+  emin = mpfr_get_emin ();
+  emax = mpfr_get_emax ();
+  set_emin (MPFR_EMIN_MIN);
+  set_emax (MPFR_EMAX_MAX);
+  if (mpfr_get_emin () != emin || mpfr_get_emax () != emax)
+    {
+      ext = 1;
+      test_small ();
+      test_large_small ();
+      check_overflow ();
+      set_emin (emin);
+      set_emax (emax);
     }
+}
 
-  mpfr_clear (x1);
-  mpfr_clear (x2);
-  mpfr_clear (y);
-  mpfr_clear (z);
-  mpfr_clear (t);
+int
+main (int argc, char *argv[])
+{
+  tests_start_mpfr ();
+
+  special ();
 
   test_large ();
-  test_large_small ();
+  alltst ();
+
+  test_generic (2, 100, 10);
 
   tests_end_mpfr ();
   return 0;