handled the case where real and imaginary part have very different

exponents in mpc_sqr


git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@34 211d60ee-9f03-0410-a15a-8952a2c7a4e4

1 files changed, 22 insertions, 1 deletions

diff --git a/sqr.c b/sqr.c
index 31c09c2..c104501 100644
--- a/sqr.c
+++ b/sqr.c
@@ -1,6 +1,6 @@
 /* mpc_sqr -- Square a complex number.
 
-Copyright (C) 2002 Free Software Foundation, Inc.
+Copyright (C) 2002 Andreas Enge, Paul Zimmermann
 
 This file is part of the MPC Library.
 
@@ -55,6 +55,27 @@ mpc_sqr (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd)
       inex_im = mpfr_set_ui (MPC_IM(a), 0, GMP_RNDN);
       return MPC_INEX(inex_re, inex_im);
    }
+   /* If the real and imaginary part of the argument have a very different */
+   /* exponent, it is not reasonable to use Karatsuba squaring; compute    */
+   /* exactly with the standard formulae instead, even if this means an    */
+   /* additional multiplication.                                           */
+   if (SAFE_ABS (mp_exp_t, MPFR_EXP (MPC_RE (b)) - MPFR_EXP (MPC_IM (b)))
+       > MPC_MAX_PREC (b) / 2)
+   {
+      mpfr_init2 (u, 2*MPFR_PREC (MPC_RE (b)));
+      mpfr_init2 (v, 2*MPFR_PREC (MPC_IM (b)));
+      
+      mpfr_sqr (u, MPC_RE (b), GMP_RNDN);
+      mpfr_sqr (v, MPC_IM (b), GMP_RNDN); /* both are exact */
+      inex_im = mpfr_mul (a->im, b->re, b->im, MPC_RND_IM (rnd));
+      mpfr_mul_2exp (a->im, a->im, 1, GMP_RNDN);
+      inex_re = mpfr_sub (a->re, u, v, MPC_RND_RE (rnd));
+      
+      mpfr_clear (u);
+      mpfr_clear (v);
+      return MPC_INEX (inex_re, inex_im);
+   }
+   
 
    mpfr_init (u);
    mpfr_init (v);