new, trivial fr_div

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

2 files changed, 16 insertions, 99 deletions

diff --git a/TODO b/TODO
index 1032657..9ba8a0d 100644
--- a/TODO
+++ b/TODO
@@ -1,4 +1,10 @@
+Efficiency:
+- from Andreas Enge 16 September 2008:
+  Once mpc_sin_cos exists, improve mpc_tan to use it
+
 New functions to implement:
+- from Andreas Enge 16 September 2008:
+  mpc_sin_cos; needs mpfr_sinh_cosh
 - from Mickael Gastineau <gastineau@imcce.fr> 14 Apr 2008:
   mpc_fma: d=a*b+c
 - from Andreas Enge 9 April 2008:
@@ -16,6 +22,7 @@ New functions to implement:
   2) cproj at any moment (relatively easy)
 - from Andreas Enge and Philippe Théveny 17 July 2008
   agm (and complex logarithm with agm ?)
+
 New tests to add:
 - from Andreas Enge and Philippe Théveny 9 April 2008
   systematic testing whether ok when output variable equals some
diff --git a/src/fr_div.c b/src/fr_div.c
index fc8a45e..771683b 100644
--- a/src/fr_div.c
+++ b/src/fr_div.c
@@ -1,6 +1,6 @@
 /* mpc_fr_div -- Divide a floating-point number by a complex number.
 
-Copyright (C) 2008 Philippe Th\'eveny
+Copyright (C) 2008 Andreas Enge
 
 This file is part of the MPC Library.
 
@@ -27,106 +27,16 @@ MA 02111-1307, USA. */
 int
 mpc_fr_div (mpc_ptr a, mpfr_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
 {
-  mpc_t cc;
-  mpc_t z;
-  mpfr_t d;
-  mp_prec_t prec;
-  mp_rnd_t zre_rnd, zim_rnd;
-  int inexact_prod, inexact_norm, inexact_re, inexact_im, loops = 0;
-  int ok_re = 0, ok_im = 0;
+   mpc_t bc;
+   int inexact;
 
-  mpfr_init (d);
-  mpc_init (z);
+   MPC_RE (bc)[0] = b [0];
+   mpfr_init (MPC_IM (bc));
+   mpfr_set_ui (MPC_IM (bc), 0, GMP_RNDN);
 
-  /* cc=conj(c), exact operation, allocating no new memory */
-  cc[0] = c[0];
-  MPFR_CHANGE_SIGN (MPC_IM (cc));
+   inexact = mpc_div (a, bc, c, rnd);
 
-  /* choose rounding direction so that |d| <= |b|/norm(c) */
-  /*   d_rnd = MPFR_SIGN (b) < 0 ? GMP_RNDU : GMP_RNDD; */
-  /* choose rounding direction so that |RE(z)| >= |RE(b/c)| and
-     |IM(z)| >= |IM(b/c)| */
-  zre_rnd = MPFR_SIGN (MPC_RE (cc)) == MPFR_SIGN (b) ? GMP_RNDU : GMP_RNDD;
-  zim_rnd = MPFR_SIGN (MPC_IM (cc)) == MPFR_SIGN (b) ? GMP_RNDU : GMP_RNDD;
+   mpfr_clear (MPC_IM (bc));
 
-  prec = MPC_MAX_PREC (a);
-
-  do
-    {
-      loops ++;
-      prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 5 : prec / 2;
-      
-      mpc_set_prec (z, prec);
-      mpfr_set_prec (d, prec);
-      
-      /* first compute norm(c)^2 */
-      inexact_norm = mpc_norm (d, c, GMP_RNDD);
-
-      /* now compute b*conjugate(c) */
-      /* We need rounding away from zero for both the real and the imagin-  */
-      /* ary part; then the final result is also rounded away from zero.    */
-      /* The error is less than 1 ulp. Since this is not implemented in     */
-      /* mpfr, we round towards zero and add 1 ulp to the absolute values   */
-      /* if they are not exact. */
-      inexact_prod = mpc_mul_fr (z, cc, b, MPC_RNDZZ);
-      inexact_re = MPC_INEX_RE (inexact_prod);
-      inexact_im = MPC_INEX_IM (inexact_prod);
-      if (inexact_re != 0)
-        mpfr_add_one_ulp (MPC_RE (z), GMP_RNDN);
-      if (inexact_im != 0)
-        mpfr_add_one_ulp (MPC_IM (z), GMP_RNDN);
-
-      /* divide the product by the norm */
-      if (inexact_norm == 0 && (inexact_re == 0 || inexact_im == 0))
-        {
-          /* The division has good chances to be exact in at least one part.*/
-          /* Since this can cause problems when not rounding to the nearest,*/
-          /* we use the division code of mpfr, which handles the situation. */
-          inexact_re |= mpfr_div (MPC_RE (z), MPC_RE (z), d, zre_rnd);
-          ok_re = mpfr_can_round (MPC_RE (z), prec - 4, zre_rnd, 
-                                  MPC_RND_RE (rnd), MPFR_PREC (MPC_RE (a)));
-
-          if (ok_re || !inexact_re) /* compute imaginary part */
-            {
-              inexact_im |= mpfr_div (MPC_IM (z), MPC_IM (z), d, zim_rnd);
-              ok_im = mpfr_can_round (MPC_IM (z), prec - 4, zim_rnd,
-                                      MPC_RND_IM (rnd), MPFR_PREC(MPC_IM(a)));
-            }
-        }
-      else
-        {
-          /* The division is inexact, so for efficiency reasons we invert q */
-          /* only once and multiply by the inverse.                         */
-          /* We do not decide about the sign of the difference. */
-          if (mpfr_ui_div (d, 1, d, GMP_RNDU))
-            {
-              /* if 1/d is inexact, the approximations of the real and
-                 imaginary part below will be inexact, unless RE(z)
-                 or IM(z) is zero */
-              inexact_re = inexact_re || (MPFR_IS_ZERO(MPC_RE(z)) == 0);
-              inexact_im = inexact_im || (MPFR_IS_ZERO(MPC_IM(z)) == 0);
-            }
-          inexact_re = mpfr_mul (MPC_RE (z), MPC_RE (z), d, zre_rnd)
-            || inexact_re;
-          ok_re = mpfr_can_round (MPC_RE (z), prec - 4, zre_rnd,
-                                  MPC_RND_RE(rnd), MPFR_PREC(MPC_RE(a)));
-
-          if (ok_re) /* compute imaginary part */
-            {
-              inexact_im = mpfr_mul (MPC_IM (z), MPC_IM (z), d, zim_rnd)
-                || inexact_im;
-              ok_im = mpfr_can_round (MPC_IM (z), prec - 4, zim_rnd,
-                                      MPC_RND_IM(rnd), MPFR_PREC(MPC_IM(a)));
-            }
-        }
-    }
-  while ((!ok_re && inexact_re) || (!ok_im && inexact_im));
-
-  mpc_set (a, z, rnd);
-
-  mpc_clear (z);
-  mpfr_clear (d);
-
-  return (MPC_INEX (inexact_re, inexact_im));
-  /* Only exactness vs. inexactness is tested, not the sign. */
+   return inexact;
 }