diff options
author | thevenyp <thevenyp@280ebfd0-de03-0410-8827-d642c229c3f4> | 2008-07-30 15:13:42 +0000 |
---|---|---|
committer | thevenyp <thevenyp@280ebfd0-de03-0410-8827-d642c229c3f4> | 2008-07-30 15:13:42 +0000 |
commit | 4452d56c95fc6c5878513b6e939f2e00cb035107 (patch) | |
tree | d2f61bb2d52d290940301d3089b253e406d6f938 | |
parent | b15ef22b767ffcbb1bc4aa50e8fbee70bf6d5679 (diff) | |
download | mpfr-4452d56c95fc6c5878513b6e939f2e00cb035107.tar.gz |
algorithms.tex: Prove the correctness of the algorithm used for mpfr_hypot
when the difference of inputs' exponents is less then exp_max - 2.
hypot.c: Change algorithm according to its description in algorithms.tex
tests/thypot.c: Fix tests (some were present but didn't trigger any error).
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@5464 280ebfd0-de03-0410-8827-d642c229c3f4
-rw-r--r-- | algorithms.tex | 202 | ||||
-rw-r--r-- | hypot.c | 75 | ||||
-rw-r--r-- | tests/thypot.c | 27 |
3 files changed, 155 insertions, 149 deletions
diff --git a/algorithms.tex b/algorithms.tex index e88cf912f..8c16e71d0 100644 --- a/algorithms.tex +++ b/algorithms.tex @@ -2267,117 +2267,119 @@ The \texttt{mpfr\_hypot} function implements the euclidean distance function: \textnormal{hypot} (x,y) = \sqrt{x^2+y^2}. \] If one of the variable is zero, then hypot is computed using the absolute -value of the other variable. -Assume that $|x| \geq |y| > 0$. Using the first degree Taylor polynomial, we -have +value of the other variable. Assume that $0 < y \leq x$. Using the first +degree Taylor polynomial, we have \[ -0 \leq \sqrt{x^2+y^2}-|x| \leq \frac{y^2}{2|x|}. +0 < \sqrt{x^2+y^2}-x < \frac{y^2}{2x}. \] -Let $p_x$ the precision of input variable $x$, $p_z$ the output precision and -$z=\circ_{p_z}(\sqrt{x^2+y^2})$ the expected result. -When $\Exp(x) - \Exp(y) > \frac{p_z+1}{2}$, we have -$\frac{y^2}{2|x|} < \frac{1}{2}\ulp_{p_z}(x)$. -When $p_z < p_x$, the stronger condition $\Exp(x) - \Exp(y) > \frac{p_x+1}{2}$ -ensures that $\sqrt{x^2+y^2} - |x| < \frac{1}{2} \ulp(x)$. -In both cases, the inequalities show that $z=\circ_{p_z}(|x|)$. +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 the 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 $x < \textnormal{hypot}(x,y)$. + +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: -When none of the above conditions is satisfied, we use the following -algorithm: \begin{quote} - Algorithm {\tt hypot}\\ - Input: $x$ and $y$ with $|y| \leq |x|$, $p$ the working precision\\ - Output: $\sqrt{x^2+y^2}$ with $p-2$ bits of precision\\ - \begin{tabular}{l p{5em} l c c l} - $u \leftarrow \Z(x^2)$ & & $\error(u)$ & $\leq$ & 1 & $\ulp(u)$\\ - $v \leftarrow \Z(y^2)$ & & $\error(v)$ & $\leq$ & 1 & $\ulp(v)$\\ - $w \leftarrow \Z(u+v)$ & & $\error(w)$ & $\leq$ & 3 & $\ulp(w)$\\ - $t \leftarrow \Z(\sqrt{w})$ & & $\error(t)$ & $\leq$ & 4 & $\ulp(t)$\\ - \end{tabular} + 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 the above algorithm, $\Z(x^2)$ or $\Z(y^2)$ may overflow though the result -$z$ is representable. Let us assume, as it is the case in MPFR, that the -minimal and maximal acceptable exponents (respectively $e_{min}$ and -$e_{max}$) verify $e_{max} = -e_{min} > 2$ and that the maximal precision is -less than $2e_{max}$. Let $s=\Z(\frac{\Exp(x) + \Exp(y) + 1}{2})$ (integer -division with rounding towards zero), we have $\frac{\Exp(x)+\Exp(y)}{2} \leq -s \leq \frac{\Exp(x) + \Exp(y)}{2} + 1$. In order to avoid overflow, we shift -input's exponents by $-s$ before entering the algorithm and shift back the -output's exponent by $+s$. Let us see now that no overflow goes on when the -exact value $\sqrt{x^2+y^2}$ is representable. - - -Let us prove first that the exponent shifts do no cause overflow. Note that -because $|y| \leq |x|$, we have $\Exp(y) \leq \Exp(x)$, thus $\Exp(y) -s \leq -\frac{\Exp(y) - \Exp(x)}{2} \leq 0$ and we just need to prove that -$\Exp(x) - s \leq e_{max}$. - -The case $\Exp(x)=-\Exp(y)$ is straightforward: $s=0$. - -Assume now that $-\Exp(y) < \Exp(x)$, then using the definition of $s$ we have -$s \geq \frac{\Exp(x) + \Exp(y)}{2} > 0$, this gives $\Exp(x) - s < \Exp(x)$. -As $e_{max}$ is the upper bound for exponents, we have $\Exp(x) -s < e_{max}$. - -In the case $\Exp(x) < -\Exp(y)$, the definition of $s$ implies that -$\frac{\Exp(x) + \Exp(y)}{2} \leq s$, so $\Exp(x)-s \leq -\frac{\Exp(x) - \Exp(y)}{2}$. But in the present case, we have -$\frac{\Exp(x) - \Exp(y)}{2} < -\Exp(y)$ and, using the lower bound for -exponents $-e_{max} \leq \Exp(y)$, we can conclude that -$\Exp(x) - s < e_{max}$. - - -Secondly, we need to show that the squares of the shifted inputs do never -overflow. Let us point out that $\Exp(x) - \Exp(y) \leq \frac{p_z + 1}{2}$, -otherwise, that case would that case would satisfy the conditions given above -and the result $z$ would have been computed using $|x|$. -Using the condition on maximal precision $p_z < 2e_{max}$, we can write -$\frac{p_z + 1}{2} \leq e_{max}$ because if $p_z$ is even, -$p_z + 1 < 2e_{max}$ and if $p_z$ is odd, $p_z +1 \leq 2e_{max}$. Thus, -$\Exp(x) - \Exp(y) \leq e_{max}$. Whatever the sign of $s$ we have -$\Exp(x) - s \leq \frac{\Exp(x) - \Exp(y)}{2}$ which is bounded by -$\frac{e_{max}}{2}$. Let $x_s = x/2^s$ and $y_s = y/2^s$ be the shifted -inputs, from the preceding inequalities we can see that $2\Exp(x_s) \leq -e_{max}$ and we can conclude that $\Exp(y_s^2) \leq \Exp(x_s^2) \leq e_{max}$, -that is to say the squares do not overflow. - - -Thirdly, noticing that the square of $y_s$ may underflow, let us prove that -the computed value is still correct despite this underflow. Actually, -$\frac{\Exp(x)+\Exp(y)}{2} \leq s \leq \frac{\Exp(x) + \Exp(y)}{2} + 1$ shows -that $-\Exp(x_s) \geq \Exp(y_s) \geq -\Exp(x_s) - 2$. [TODO] - - -Fourthly, let us see that the addition does not overflow. The addition -overflows when $\Exp(w) > e_{max}$, since $u$ and $v$ are representable and -the rounding mode is towards zero, the only case in which this happens is when -$u + v \geq 2^{\Exp(u)}$ and $\Exp(w) = \Exp(u) + 1$. -Then $v \geq 2^{\Exp(u)} - u$ and $\Exp(u) = e_{max}$. -We proved above that $2\Exp(x_s) \leq e_{max}$ and we know that -$\Exp(u) \leq 2\Exp(x_s)$, so using the last three relations we can write -$2\Exp(x_s) = \Exp(u) = e_{max}$ which is not possible if $e_{max}$ is odd. -Assume that $e_{max}$ is even, then $u \leq (1-2^{-p})2^{\Exp(u)}$ where $p$ -is the working precision in the algorithm and we have $v \geq 2^{\Exp(u)-p}$, -on the other hand, because of the rounding towards zero mode -$2^{2\Exp(y_s)} \geq v$, thus $p \geq 2\left(\Exp(x_s) - \Exp(y_s)\right)$. -Using the fact that $\Exp(y_s) \geq -\Exp(x_s) - 2$, we can write $p \geq -4\Exp(x_s) + 4 > 2e_{max}$ which is incompatible with the upper limit for -precision. In conclusion, the addition does not overflow in any case. +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 we have +seen 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 not representable. -As $\sqrt{x^2+y^2} \geq |x|$, \texttt{mpfr\_hypot} should never underflow. -We have $\Exp(x_s) = \Exp(x) - s$ by definition of $s$, -$2\Exp(x) - 2s \geq \Exp(u) \geq 2\Exp(x) - 2s - 1$ because the squared -mantissa of $x$ may be less than $\frac{1}{2}$, -$2\Exp(x) - 2s + 1 \geq \Exp(w) \geq 2\Exp(x) - 2s - 1$ because $u \geq v$, -$\Exp(x) - s + 1 \geq \Exp(t) \geq \Exp(x) - s$ ($\Exp(\sqrt{r}) = n$ when -$\Exp(r) = 2n$ is even and $\Exp(\sqrt{r}) = n + 1$ when $\Exp(r) = 2n + 1$ is -odd), and $\Exp(x) + 1 \geq \Exp(z) \geq \Exp(x)$ which shows that it does not -unverflow but may overflow when $\Exp(x) = e_{max}$. -[TODO] +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(u)$\\ + $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, we have lost 2 bits of precision when $\max(p_x, p_y) \leq p$ and 4 bits +when $p$ does not verify this relation. \subsection{The floating multiply-add function} @@ -31,10 +31,12 @@ 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); @@ -72,26 +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). - FIXME: ^^^^^^^^^^^^^^^^^^^^^^^^^ is meaningless. - 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) + 1) / 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)) @@ -104,7 +96,7 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) } else /* GMP_RNDZ, GMP_RNDD, GMP_RNDN */ { - if (MPFR_LIKELY (Nz >= Nx)) + if (MPFR_LIKELY (Nz >= N)) { mpfr_abs (z, x, rnd_mode); /* exact */ MPFR_RET (-1); @@ -113,7 +105,7 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) { 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); @@ -128,43 +120,40 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) /* General case */ - Ny = MPFR_PREC(y); /* Precision of input variable */ + /* FIXME: the following algorithm is correct only for + diff_exp <= MPFR_EMAX_DEFAULT - 2 */ - /* 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; + N = MAX (MPFR_PREC (x), MPFR_PREC (y)); + + /* 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); - /* Scale x and y to avoid overflow/underflow in x^2 and y^2. - After scaling, we need to have exponent values as small as - possible in absolute value, for both x and y. So, the best - choice is to scale by about (Ex + Ey) / 2. */ - sh = (Ex + Ey) / 2; + /* 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_2si (te, x, sh, GMP_RNDZ); /* exact since Nt >= Nx */ - mpfr_div_2si (ti, y, sh, GMP_RNDZ); /* exact since Nt >= Ny */ - exact = mpfr_sqr (te, te, GMP_RNDZ); /* x^2 */ - exact |= mpfr_sqr (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; - /* update the precision */ MPFR_ZIV_NEXT (loop, Nt); mpfr_set_prec (t, Nt); mpfr_set_prec (te, Nt); @@ -172,7 +161,7 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) } MPFR_ZIV_FREE (loop); - MPFR_BLOCK (flags, inexact = mpfr_mul_2si (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); diff --git a/tests/thypot.c b/tests/thypot.c index de691f5ad..47507c3ed 100644 --- a/tests/thypot.c +++ b/tests/thypot.c @@ -92,19 +92,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); } |