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+ssh://scm.gforge.inria.fr/svnroot/mpc/trunk@1202 211d60ee-9f03-0410-a15a-8952a2c7a4e4

4 files changed, 112 insertions, 28 deletions

diff --git a/doc/algorithms.tex b/doc/algorithms.tex
index 790eca1..edad662 100644
--- a/doc/algorithms.tex
+++ b/doc/algorithms.tex
@@ -3,6 +3,7 @@
 \usepackage[a4paper]{geometry}
 \usepackage[utf8]{inputenc}
 \usepackage[T1]{fontenc}
+\usepackage{ae}
 \usepackage{amsmath,amssymb}
 \usepackage{hyperref}
 \usepackage{comment}
@@ -46,7 +47,7 @@
 
 \title {MPC: Algorithms and Error Analysis}
 \author {Andreas Enge \and Philippe Th\'eveny \and Paul Zimmermann}
-\date {Draft; June 27, 2012}
+\date {Draft; June 28, 2012}
 
 \begin {document}
 \maketitle
@@ -454,6 +455,37 @@ For the sine function, a completely analogous argument shows that
 \eqref {eq:proprealcos} also holds.
 
 
+\subsubsection {Logarithm}
+\label {sssec:propreallog}
+
+Let
+\[
+\appro x = \log (1 + \appro {x_1})
+\]
+for $\appro {x_1} > -1$.
+By the mean value theorem, there is a $\xi$ between $x_1$ and $\appro {x_1}$
+such that
+\[
+\error (\appro x) = \frac {1}{1 + \xi} \error (\appro {x_1})
+\leq \frac {1}{1 + \min (x_1, \appro {x_1})} \error (\appro {x_1}).
+\]
+For $x_1 > 0$, this implies
+\begin {eqnarray*}
+\error (\appro x)
+& \leq & \error (\appro {x_1})
+\leq
+k \, 2^{\Exp (\appro {x_1}) - \Exp (\appro x)}
+\, 2^{\Exp (\appro x) - p} \\
+& \leq & 2 \, k \, \frac {\appro {x_1}}{\appro x} \, 2^{\Exp (\appro x) - p} \\
+& \leq & 2 \, k \, \frac {\appro {x_1}}{\appro {x_1} - \appro {x_1}^2/2}
+\, 2^{\Exp (\appro x) - p}
+\end {eqnarray*}
+using $\log (1 + z) \geq z - z^2/2$ for $z > 0$.
+For $0 < x_1 \leq 1$, we have $\appro {x_1}^2/2 \leq \appro {x_1}/2$ and
+\[
+\error (\appro x)
+\leq 4 \, k \, 2^{\Exp (\appro x) - p}.
+\]
 
 
 \subsection {Complex functions}
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);
 }
diff --git a/tests/log.dat b/tests/log.dat
index 0f590d0..6972cfc 100644
--- a/tests/log.dat
+++ b/tests/log.dat
@@ -166,10 +166,10 @@
 + + 53 0x162E4301C0025p-23   2  0x1p0         2 0x1p67108864 2 0x1p67108864 U U
 + + 53 0xB1721802E8F1p-21    2  0x1p0       2 0x1p134217728 2 0x1p134217728 U U
 + + 53 0x2C5C85FF57581p-22   2  0x1p0       2 0x1p268435456 2 0x1p268435456 U U
-# Due to intermediate overflow, the following result has real part infinity
-# instead of the correct result. Since this may happen in other parts of the
-# library as well, we do not consider it a bug for the moment.
-# + + 53 0x58B90BFD4BCBFp-22   2  0x1p0       2 0x1p536870912 2 0x1p536870912 U U
+
+# Example leading to intermediate overflow in x^2+y^2
++ + 53 0x58B90BFD4BCBFp-22   2  0x1p0       2 0x1p536870912 2 0x1p536870912 U U
+
 # Due to intermediate "underflow" (see FIXME in log.c), the following example
 # loops.
 # log(-1 + i*eps)