log: alternative implementation that avoids intermediate overflows

     It is probably slower (two calls to log) and should be combined with
     the previous approach.
     Problem of "underflow" (log of number close to 1) not yet solved.


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

2 files changed, 75 insertions, 23 deletions

diff --git a/src/fma.c b/src/fma.c
index d4be40f..726ec4e 100644
--- a/src/fma.c
+++ b/src/fma.c
@@ -154,6 +154,20 @@ mpc_fma (mpc_ptr r, mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
   mpc_init3 (ab, wpre, wpim);
   for (i = 0; i < 2; ++i)
     {
+#define MPC_OUT(x)                                              \
+do {                                                            \
+  printf (#x "[%lu,%lu]=", (unsigned long int) MPC_PREC_RE (x), \
+      (unsigned long int) MPC_PREC_IM (x));                     \
+  mpc_out_str (stdout, 2, 0, x, MPC_RNDNN);                     \
+  printf ("\n");                                                \
+} while (0)
+
+#define MPFR_OUT(x)                                             \
+do {                                                            \
+  printf (#x "[%lu]=", (unsigned long int) mpfr_get_prec (x));  \
+  mpfr_out_str (stdout, 2, 0, x, GMP_RNDN);                     \
+  printf ("\n");                                                \
+} while (0)
       mpc_mul (ab, a, b, MPC_RNDZZ);
       if (mpfr_zero_p (mpc_realref(ab)) || mpfr_zero_p (mpc_imagref(ab)))
         break;
diff --git a/src/log.c b/src/log.c
index 5a91fb4..3b9f799 100644
--- a/src/log.c
+++ b/src/log.c
@@ -24,11 +24,14 @@ along with this program. If not, see http://www.gnu.org/licenses/ .
 int
 mpc_log (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd){
    int ok=0;
-   mpfr_t w;
+   mpfr_srcptr x, y;
+   mpfr_t v, w;
    mpfr_prec_t prec;
    int loops = 0;
    int re_cmp, im_cmp;
    int inex_re, inex_im, inex;
+   int err;
+   mpfr_exp_t expw;
 
    /* special values: NaN and infinities */
    if (!mpc_fin_p (op)) {
@@ -112,40 +115,75 @@ mpc_log (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd){
    }
 
    prec = MPC_PREC_RE(rop);
+   mpfr_init2 (v, prec);
    mpfr_init2 (w, prec);
    /* let op = x + iy; log = 1/2 log (x^2 + y^2) + i atan2 (y, x) */
-   /* loop for the real part: log (x^2 + y^2)                    */
+   /* loop for the real part:
+      compute 1/2 log (x^2 + y^2) = log |x| + 1/2 * log (1 + (y/x)^2)
+         if |x| >= |y|; otherwise, exchange x and y               */
+   if (mpfr_cmpabs (mpc_realref (op), mpc_imagref (op)) >= 0) {
+      x = mpc_realref (op);
+      y = mpc_imagref (op);
+   }
+   else {
+      x = mpc_imagref (op);
+      y = mpc_realref (op);
+   }
    do {
+#define MPC_OUT(x)                                              \
+do {                                                            \
+  printf (#x "[%lu,%lu]=", (unsigned long int) MPC_PREC_RE (x), \
+      (unsigned long int) MPC_PREC_IM (x));                     \
+  mpc_out_str (stdout, 2, 0, x, MPC_RNDNN);                     \
+  printf ("\n");                                                \
+} while (0)
+
+#define MPFR_OUT(x)                                             \
+do {                                                            \
+  printf (#x "[%lu]=", (unsigned long int) mpfr_get_prec (x));  \
+  mpfr_out_str (stdout, 2, 0, x, GMP_RNDN);                     \
+  printf ("\t");                                                 \
+  mpfr_out_str (stdout, 10, 0, x, GMP_RNDN);                     \
+  printf ("\n");                                                \
+} while (0)
       loops ++;
-      prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2;
+      prec += mpc_ceil_log2 (prec) + 4;
+      mpfr_set_prec (v, prec);
       mpfr_set_prec (w, prec);
 
-      /* w is rounded down */
-      inex = mpc_norm (w, op, GMP_RNDD);
-         /* error 1 ulp */
-      MPC_ASSERT (!mpfr_inf_p (w));
-         /* FIXME: intermediate overflow; the logarithm may be representable */
-      MPC_ASSERT (!(inex && mpfr_cmp_ui (w, 1) == 0));
-         /* FIXME: this is like an underflow; the following call to log will
-            compute 0 instead of a positive result                           */
-      
-      mpfr_log (w, w, GMP_RNDD);
-      /* generic error of log: (2^(2 - exp(w)) + 1) ulp */
+      mpfr_div (v, y, x, GMP_RNDN); /* error 0.5 ulp */
+      mpfr_sqr (v, v, GMP_RNDN);
+         /* generic error of multiplication:
+            0.5 + 2*0.5*(2+0.5*2^(1-prec)) <= 2.51... since prec >= 6 */
+      mpfr_log1p (v, v, GMP_RNDN);
+         /* error 0.5 + 4*2.51... = 10.54... , see algorithms.tex */
+      mpfr_div_2ui (v, v, 1, GMP_RNDN);
 
-      if (mpfr_get_exp (w) >= 2)
-         ok = mpfr_can_round (w, prec - 2, GMP_RNDD,
-            MPC_RND_RE(rnd), MPC_PREC_RE(rop));
-      else
-         ok = mpfr_can_round (w, prec - 3 + mpfr_get_exp (w), GMP_RNDD,
-            MPC_RND_RE(rnd), MPC_PREC_RE(rop));
-   } while (ok == 0);
+      mpfr_abs (w, x, GMP_RNDN); /* exact */
+      mpfr_log (w, w, GMP_RNDN); /* error 0.5 ulp */
+      expw = mpfr_get_exp (w);
 
+      mpfr_add (w, w, v, GMP_RNDN);
+      if (!mpfr_signbit (w)) /* v is positive, so no cancellation;
+                                error 11.54... ulp; error counts lost bits */
+         err = 4;
+      else
+         err =   MPC_MAX (4 + mpfr_get_exp (v),
+                 /* 10.54... ulp (v) rewritten in ulp (result, now in w) */
+                          -1 + expw             - mpfr_get_exp (w)
+                 /* 0.5 ulp (previous w), rewritten in ulp (result) */
+                 ) + 2;
+   
+   } while (!mpfr_can_round (w, prec - err, GMP_RNDN, GMP_RNDZ,
+            mpfr_get_prec (mpc_realref (rop)) + (MPC_RND_RE (rnd) == GMP_RNDN)));
+      
    /* imaginary part */
    inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), mpc_realref (op),
                          MPC_RND_IM (rnd));
 
-   /* set the real part; cannot be done before when rop==op */
-   inex_re = mpfr_div_2ui (mpc_realref(rop), w, 1ul, MPC_RND_RE (rnd));
+   /* set the real part; cannot be done before if rop==op */
+   inex_re = mpfr_set (mpc_realref (rop), w, MPC_RND_RE (rnd));
+   mpfr_clear (v);
    mpfr_clear (w);
    return MPC_INEX(inex_re, inex_im);
 }