From bef243d0bdd8914b1f7a364b1dd32f3df882c389 Mon Sep 17 00:00:00 2001 From: enge Date: Thu, 1 Mar 2012 21:07:48 +0000 Subject: sqr: copied code for imaginary part from naive computation also to the Karatsuba algorithm sqr.dat: added examples with imaginary part previously miscomputed by Karatsuba git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@1133 211d60ee-9f03-0410-a15a-8952a2c7a4e4 --- src/sqr.c | 32 +++++++++++--------------------- tests/sqr.dat | 4 ++++ 2 files changed, 15 insertions(+), 21 deletions(-) diff --git a/src/sqr.c b/src/sqr.c index 109e314..47a3f7e 100644 --- a/src/sqr.c +++ b/src/sqr.c @@ -180,6 +180,7 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_prec_t prec; int inex_re, inex_im, inexact; mpfr_exp_t emin, emax; + int saved_underflow; /* special values: NaN and infinities */ if (!mpc_fin_p (op)) { @@ -252,20 +253,10 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) exactly with the standard formulae instead, even if this means an additional multiplication. Using the approach copied from mul, over- and underflows are also handled correctly. */ - int saved_underflow; inex_re = mpfr_fmma (rop->re, x, x, op->im, op->im, -1, MPC_RND_RE (rnd)); - saved_underflow = mpfr_underflow_p (); - mpfr_clear_underflow (); - inex_im = mpfr_mul (rop->im, x, op->im, MPC_RND_IM (rnd)); - if (!mpfr_underflow_p ()) - inex_im |= mpfr_mul_2exp (rop->im, rop->im, 1, MPC_RND_IM (rnd)); - /* We must not multiply by 2 if rop->im has been set to the smallest - representable number. */ - if (saved_underflow) - mpfr_set_underflow (); } else { /* Karatsuba squaring */ @@ -367,19 +358,18 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_clear (u); mpfr_clear (v); - - if (mpfr_get_exp (x) + mpfr_get_exp(mpc_imagref (op)) <= emin + 1) { - mpfr_mul_2ui (x, x, 1, GMP_RNDN); - inex_im = mpfr_mul (mpc_imagref (rop), x, mpc_imagref (op), MPC_RND_IM (rnd)); - } - else { - inex_im = mpfr_mul (mpc_imagref (rop), x, mpc_imagref (op), MPC_RND_IM (rnd)); - /* checks that no underflow occurs (an overflow can occur here) */ - MPC_ASSERT (mpfr_zero_p (mpc_imagref (rop)) == 0); - mpfr_mul_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN); - } } + saved_underflow = mpfr_underflow_p (); + mpfr_clear_underflow (); + inex_im = mpfr_mul (rop->im, x, op->im, MPC_RND_IM (rnd)); + if (!mpfr_underflow_p ()) + inex_im |= mpfr_mul_2exp (rop->im, rop->im, 1, MPC_RND_IM (rnd)); + /* We must not multiply by 2 if rop->im has been set to the smallest + representable number. */ + if (saved_underflow) + mpfr_set_underflow (); + if (rop == op) mpfr_clear (x); diff --git a/tests/sqr.dat b/tests/sqr.dat index ee59b0e..9d56f0e 100644 --- a/tests/sqr.dat +++ b/tests/sqr.dat @@ -160,6 +160,10 @@ # another interesting one with not exactly the same behaviour - - 100 -inf 100 -inf 100 -0xf@192058806 100 0x1@192058873 N N 0 + 100 0 100 inf 100 0x1@192058806 100 0x1@192058806 N N +# Re(op)*Im(op) can be computed, but multiplication by 2 triggers overflow +0 + 100 0 100 inf 100 0b1@536870911 100 0b1@536870911 N N +0 - 10 0 10 0b1.111111111e1073741822 100 0b1@536870911 100 0b1@536870911 N D +0 - 10 0 10 0b1.111111111e1073741822 100 0b1@536870912 100 0b1@536870912 N D 0 0 10 0 10 0b1e-1073741823 100 0b1@-536870912 100 0b1@-536870912 N N 0 - 10 0 10 0 100 0b1@-536870913 100 0b1@-536870913 N N 0 + 10 0 10 0b1@-1073741824 100 0b1@-536870913 100 0b1@-536870913 N U -- cgit v1.2.1