added new function mpc_log10

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

9 files changed, 468 insertions, 7 deletions

diff --git a/NEWS b/NEWS
index 980f19e..0b3887a 100644
--- a/NEWS
+++ b/NEWS
@@ -4,6 +4,7 @@ Changes in version 1.0:
   - 0^0, which returned (NaN,NaN) previously, now returns (1,+0).
   - removed compatibility with K&R compilers, untestable due to lack of
     such compilers
+  - New function mpc_log10.
   - Bug fixes:
     - div and norm now return a value indicating the effective rounding
       direction, as the other functions.
diff --git a/doc/mpc.texi b/doc/mpc.texi
index acb4576..a724877 100644
--- a/doc/mpc.texi
+++ b/doc/mpc.texi
@@ -959,11 +959,12 @@ rounded according to @var{rnd} with the precision of @var{rop}.
 @end deftypefun
 
 @deftypefun int mpc_log (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd})
-Set @var{rop} to the logarithm of @var{op},
+@deftypefunx int mpc_log10 (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd})
+Set @var{rop} to the natural and base-10 logarithm of @var{op} respectively,
 rounded according to @var{rnd} with the precision of @var{rop}.
 The principal branch is chosen, with the branch cut on the negative real axis,
 so that the imaginary part of the result lies in
-@math{]-\pi , \pi]}.
+@math{]-\pi , \pi]} and @math{]-\pi/log(10) , \pi/log(10)]} respectively.
 @end deftypefun
 
 
diff --git a/src/Makefile.am b/src/Makefile.am
index 0f0ddd6..5cea5c2 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -23,7 +23,7 @@ libmpc_la_SOURCES = mpc-impl.h abs.c acos.c acosh.c add.c add_fr.c         \
   add_si.c add_ui.c arg.c asin.c asinh.c atan.c atanh.c clear.c cmp.c      \
   cmp_si_si.c conj.c cos.c cosh.c div_2exp.c div.c div_fr.c div_ui.c exp.c \
   fma.c fr_div.c fr_sub.c get_prec2.c get_prec.c get_version.c get_x.c     \
-  imag.c init2.c init3.c inp_str.c log.c mem.c mul_2exp.c mul.c            \
+  imag.c init2.c init3.c inp_str.c log.c log10.c mem.c mul_2exp.c mul.c    \
   mul_fr.c mul_i.c mul_si.c mul_ui.c neg.c norm.c out_str.c pow.c pow_fr.c \
   pow_ld.c pow_d.c pow_si.c pow_ui.c pow_z.c proj.c real.c urandom.c set.c \
   set_prec.c set_str.c set_x.c set_x_x.c sin.c sin_cos.c sinh.c sqr.c      \
diff --git a/src/log10.c b/src/log10.c
new file mode 100644
index 0000000..88dafd6
--- /dev/null
+++ b/src/log10.c
@@ -0,0 +1,248 @@
+/* mpc_log10 -- Take the base-10 logarithm of a complex number.
+
+Copyright (C) 2012 INRIA
+
+This file is part of GNU MPC.
+
+GNU MPC is free software; you can redistribute it and/or modify it under
+the terms of the GNU Lesser General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
+more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program. If not, see http://www.gnu.org/licenses/ .
+*/
+
+#include "mpc-impl.h"
+
+/* Auxiliary functions which implement Ziv's strategy for special cases.
+   if flag = 0: compute only real part
+   if flag = 1: compute only imaginary
+   Exact cases should be dealt with separately. */
+static int
+mpc_log10_aux (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd, int flag, int nb)
+{
+  mp_prec_t prec = (MPFR_PREC_MIN > 4) ? MPFR_PREC_MIN : 4;
+  mpc_t tmp;
+  mpfr_t log10;
+  int ok = 0, ret;
+
+  prec = mpfr_get_prec ((flag == 0) ? mpc_realref (rop) : mpc_imagref (rop));
+  prec += 10;
+  mpc_init2 (tmp, prec);
+  mpfr_init2 (log10, prec);
+  while (ok == 0)
+    {
+      mpfr_set_ui (log10, 10, GMP_RNDN); /* exact since prec >= 4 */
+      mpfr_log (log10, log10, GMP_RNDN);
+      /* In each case we have two roundings, thus the final value is
+         x * (1+u)^2 where x is the exact value, and |u| <= 2^(-prec-1).
+         Thus the error is always less than 3 ulps. */
+      switch (nb)
+        {
+        case 0: /* imag <- atan2(y/x) */
+          mpfr_atan2 (mpc_imagref (tmp), mpc_imagref (op), mpc_realref (op),
+                      MPC_RND_IM (rnd));
+          mpfr_div (mpc_imagref (tmp), mpc_imagref (tmp), log10, GMP_RNDN);
+          ok = mpfr_can_round (mpc_imagref (tmp), prec - 2, GMP_RNDZ,
+                               GMP_RNDZ, MPC_PREC_IM(rop) +
+                               MPC_RND_IM (rnd) == GMP_RNDN);
+          if (ok)
+            ret = mpfr_set (mpc_imagref (rop), mpc_imagref (tmp),
+                            MPC_RND_IM (rnd));
+          break;
+        case 1: /* real <- log(x) */
+          mpfr_log (mpc_realref (tmp), mpc_realref (op), MPC_RND_RE (rnd));
+          mpfr_div (mpc_realref (tmp), mpc_realref (tmp), log10, GMP_RNDN);
+          ok = mpfr_can_round (mpc_realref (tmp), prec - 2, GMP_RNDZ,
+                               GMP_RNDZ, MPC_PREC_RE(rop) +
+                               MPC_RND_RE (rnd) == GMP_RNDN);
+          if (ok)
+            ret = mpfr_set (mpc_realref (rop), mpc_realref (tmp),
+                            MPC_RND_RE (rnd));
+          break;
+        case 2: /* imag <- pi */
+          mpfr_const_pi (mpc_imagref (tmp), MPC_RND_IM (rnd));
+          mpfr_div (mpc_imagref (tmp), mpc_imagref (tmp), log10, GMP_RNDN);
+          ok = mpfr_can_round (mpc_imagref (tmp), prec - 2, GMP_RNDZ,
+                               GMP_RNDZ, MPC_PREC_IM(rop) +
+                               MPC_RND_IM (rnd) == GMP_RNDN);
+          if (ok)
+            ret = mpfr_set (mpc_imagref (rop), mpc_imagref (tmp),
+                            MPC_RND_IM (rnd));
+          break;
+        case 3: /* real <- log(y) */
+          mpfr_log (mpc_realref (tmp), mpc_imagref (op), MPC_RND_RE (rnd));
+          mpfr_div (mpc_realref (tmp), mpc_realref (tmp), log10, GMP_RNDN);
+          ok = mpfr_can_round (mpc_realref (tmp), prec - 2, GMP_RNDZ,
+                               GMP_RNDZ, MPC_PREC_RE(rop) +
+                               MPC_RND_RE (rnd) == GMP_RNDN);
+          if (ok)
+            ret = mpfr_set (mpc_realref (rop), mpc_realref (tmp),
+                            MPC_RND_RE (rnd));
+          break;
+        }
+      prec += prec / 2;
+      mpc_set_prec (tmp, prec);
+      mpfr_set_prec (log10, prec);
+    }
+  mpc_clear (tmp);
+  mpfr_clear (log10);
+  return ret;
+}
+
+int
+mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
+{
+  int ok = 0, loops = 0, re_cmp, im_cmp, inex_re, inex_im, negative_zero;
+  mpfr_t w;
+  mpfr_prec_t prec;
+  mpfr_rnd_t rnd_im;
+  mpc_t ww;
+  mpc_rnd_t invrnd;
+
+  /* special values: NaN and infinities: same as mpc_log */
+  if (!mpc_fin_p (op)) /* real or imaginary parts are NaN or Inf */
+    {
+      if (mpfr_nan_p (mpc_realref (op)))
+        {
+          if (mpfr_inf_p (mpc_imagref (op)))
+            /* (NaN, Inf) -> (+Inf, NaN) */
+            mpfr_set_inf (mpc_realref (rop), +1);
+          else
+            /* (NaN, xxx) -> (NaN, NaN) */
+            mpfr_set_nan (mpc_realref (rop));
+         mpfr_set_nan (mpc_imagref (rop));
+         inex_im = 0; /* Inf/NaN is exact */
+        }
+      else if (mpfr_nan_p (mpc_imagref (op)))
+        {
+          if (mpfr_inf_p (mpc_realref (op)))
+            /* (Inf, NaN) -> (+Inf, NaN) */
+            mpfr_set_inf (mpc_realref (rop), +1);
+          else
+            /* (xxx, NaN) -> (NaN, NaN) */
+            mpfr_set_nan (mpc_realref (rop));
+          mpfr_set_nan (mpc_imagref (rop));
+          inex_im = 0; /* Inf/NaN is exact */
+        }
+      else /* We have an infinity in at least one part. */
+        {
+          /* (+Inf, y) -> (+Inf, 0) for finite positive-signed y */
+          if (mpfr_inf_p (mpc_realref (op)) && mpfr_signbit (mpc_realref (op))
+              == 0 && mpfr_number_p (mpc_imagref (op)))
+            inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op),
+                                  mpc_realref (op), MPC_RND_IM (rnd));
+          else
+            /* (xxx, Inf) -> (+Inf, atan2(Inf/xxx))
+               (Inf, yyy) -> (+Inf, atan2(yyy/Inf)) */
+            inex_im = mpc_log10_aux (rop, op, rnd, 1, 0);
+          mpfr_set_inf (mpc_realref (rop), +1);
+        }
+      return MPC_INEX(0, inex_im);
+    }
+
+   /* special cases: real and purely imaginary numbers */
+  re_cmp = mpfr_cmp_ui (mpc_realref (op), 0);
+  im_cmp = mpfr_cmp_ui (mpc_imagref (op), 0);
+  if (im_cmp == 0) /* Im(op) = 0 */
+    {
+      if (re_cmp == 0) /* Re(op) = 0 */
+        {
+          if (mpfr_signbit (mpc_realref (op)) == 0)
+            inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op),
+                                  mpc_realref (op), MPC_RND_IM (rnd));
+          else
+            inex_im = mpc_log10_aux (rop, op, rnd, 1, 0);
+          mpfr_set_inf (mpc_realref (rop), -1);
+          inex_re = 0; /* -Inf is exact */
+        }
+      else if (re_cmp > 0)
+        {
+          inex_re = mpfr_log10 (mpc_realref (rop), mpc_realref (op),
+                                MPC_RND_RE (rnd));
+          inex_im = mpfr_set (mpc_imagref (rop), mpc_imagref (op),
+                              MPC_RND_IM (rnd));
+        }
+      else /* log10(x + 0*i) for negative x */
+        { /* op = x + 0*i; let w = -x = |x| */
+          negative_zero = mpfr_signbit (mpc_imagref (op));
+          if (negative_zero)
+            rnd_im = INV_RND (MPC_RND_IM (rnd));
+          else
+            rnd_im = MPC_RND_IM (rnd);
+          ww->re[0] = *mpc_realref (op);
+          MPFR_CHANGE_SIGN (ww->re);
+          ww->im[0] = *mpc_imagref (op);
+          if (mpfr_cmp_ui (ww->re, 1) == 0)
+            inex_re = mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd));
+          else
+            inex_re = mpc_log10_aux (rop, ww, rnd, 0, 1);
+          inex_im = mpc_log10_aux (rop, op, RNDC(0,rnd_im), 1, 2);
+          if (negative_zero)
+            {
+              mpc_conj (rop, rop, MPC_RNDNN);
+              inex_im = -inex_im;
+            }
+        }
+      return MPC_INEX(inex_re, inex_im);
+    }
+   else if (re_cmp == 0)
+     {
+       if (im_cmp > 0)
+         {
+           inex_re = mpc_log10_aux (rop, op, rnd, 0, 3);
+           inex_im = mpc_log10_aux (rop, op, rnd, 1, 2);
+           /* division by 2 does not change the ternary flag */
+           mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN);
+         }
+       else
+         {
+           ww->re[0] = *mpc_realref (op);
+           ww->im[0] = *mpc_imagref (op);
+           MPFR_CHANGE_SIGN (ww->im);
+           inex_re = mpc_log10_aux (rop, ww, rnd, 0, 3);
+           invrnd = RNDC(0, INV_RND (MPC_RND_IM (rnd)));
+           inex_im = mpc_log10_aux (rop, op, invrnd, 1, 2);
+           /* division by 2 does not change the ternary flag */
+           mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN);
+           mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), GMP_RNDN);
+           inex_im = -inex_im; /* negate the ternary flag */
+         }
+       return MPC_INEX(inex_re, inex_im);
+     }
+
+  /* generic case: neither Re(op) nor Im(op) is NaN, Inf or zero */
+  prec = MPC_PREC_RE(rop);
+  mpfr_init2 (w, prec);
+  mpc_init2 (ww, prec);
+  /* let op = x + iy; compute log(op)/log(10) */
+  while (ok == 0)
+    {
+      loops ++;
+      prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2;
+      mpfr_set_prec (w, prec);
+      mpc_set_prec (ww, prec);
+
+      mpc_log (ww, op, MPC_RNDNN);
+      mpfr_set_ui (w, 10, GMP_RNDN); /* exact since prec >= 4 */
+      mpfr_log (w, w, GMP_RNDN);
+      mpc_div_fr (ww, ww, w, MPC_RNDNN);
+
+      ok = mpfr_can_round (mpc_realref (ww), prec - 2, GMP_RNDZ, GMP_RNDZ,
+                           MPC_PREC_RE(rop) + MPC_RND_RE (rnd) == GMP_RNDN);
+      ok = ok && mpfr_can_round (mpc_imagref (ww), prec-2, GMP_RNDZ, GMP_RNDZ,
+                             MPC_PREC_IM(rop) + MPC_RND_IM (rnd) == GMP_RNDN);
+    }
+
+  inex_re = mpfr_set (mpc_realref(rop), mpc_realref (ww), MPC_RND_RE (rnd));
+  inex_im = mpfr_set (mpc_imagref(rop), mpc_imagref (ww), MPC_RND_IM (rnd));
+  mpfr_clear (w);
+  mpc_clear (ww);
+  return MPC_INEX(inex_re, inex_im);
+}
diff --git a/src/mpc.h b/src/mpc.h
index ab87e7c..c98c3fb 100644
--- a/src/mpc.h
+++ b/src/mpc.h
@@ -188,6 +188,7 @@ __MPC_DECLSPEC int  mpc_cmp       (mpc_srcptr, mpc_srcptr);
 __MPC_DECLSPEC int  mpc_cmp_si_si (mpc_srcptr, long int, long int);
 __MPC_DECLSPEC int  mpc_exp       (mpc_ptr, mpc_srcptr, mpc_rnd_t);
 __MPC_DECLSPEC int  mpc_log       (mpc_ptr, mpc_srcptr, mpc_rnd_t);
+__MPC_DECLSPEC int  mpc_log10     (mpc_ptr, mpc_srcptr, mpc_rnd_t);
 __MPC_DECLSPEC int  mpc_sin       (mpc_ptr, mpc_srcptr, mpc_rnd_t);
 __MPC_DECLSPEC int  mpc_cos       (mpc_ptr, mpc_srcptr, mpc_rnd_t);
 __MPC_DECLSPEC int  mpc_sin_cos   (mpc_ptr, mpc_ptr, mpc_srcptr, mpc_rnd_t, mpc_rnd_t);
diff --git a/tests/Makefile.am b/tests/Makefile.am
index d4f425e..dede832 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -29,10 +29,10 @@ LOADLIBES=$(DEFS) -I$(top_srcdir)/src -I$(top_builddir) $(CPPFLAGS) $(CFLAGS) -L
 
 check_PROGRAMS = tabs tacos tacosh tadd tadd_fr tadd_si tadd_ui targ tasin    \
 tasinh tatan tatanh tconj tcos tcosh tdiv tdiv_2exp tdiv_fr tdiv_ui texp tfma \
-tfr_div tfr_sub timag tio_str tlog tmul tmul_2exp tmul_fr tmul_i tmul_si      \
-tmul_ui tneg tnorm tpow tpow_ld tpow_d tpow_fr tpow_si tpow_ui tpow_z tprec   \
-tproj treal treimref tset tsin tsin_cos tsinh tsqr tsqrt tstrtoc tsub tsub_fr \
-tsub_ui tswap ttan ttanh tui_div tui_ui_sub tget_version
+tfr_div tfr_sub timag tio_str tlog tlog10 tmul tmul_2exp tmul_fr tmul_i       \
+tmul_si tmul_ui tneg tnorm tpow tpow_ld tpow_d tpow_fr tpow_si tpow_ui tpow_z \
+tprec tproj treal treimref tset tsin tsin_cos tsinh tsqr tsqrt tstrtoc tsub   \
+tsub_fr tsub_ui tswap ttan ttanh tui_div tui_ui_sub tget_version
 
 check_LTLIBRARIES=libmpc-tests.la
 libmpc_tests_la_SOURCES=mpc-tests.h random.c tgeneric.c read_data.c           \
diff --git a/tests/log.dat b/tests/log.dat
index 7222fb5..73818f3 100644
--- a/tests/log.dat
+++ b/tests/log.dat
@@ -95,6 +95,7 @@
 
 # log(nan + i*inf) = +inf + i*nan
 0 0 2 +inf 2 nan                   2 nan 2 +inf  N N
+0 0 2 +inf 2 nan                   2 nan 2 -inf  N N
 
 # log(nan + i*nan) = nan + i*nan
 0 0 2 nan 2 nan                    2 nan 2 nan   N N
diff --git a/tests/log10.dat b/tests/log10.dat
new file mode 100644
index 0000000..acd24dd
--- /dev/null
+++ b/tests/log10.dat
@@ -0,0 +1,172 @@
+# Data test file for mpc_log10.
+#
+# Copyright (C) 2012 INRIA
+#
+# This file is part of GNU MPC.
+#
+# GNU MPC is free software; you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by the
+# Free Software Foundation; either version 3 of the License, or (at your
+#o ption) any later version.
+#
+# GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
+# more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see http://www.gnu.org/licenses/ .
+
+# See file sin.dat for the format description.
+
+# Special values, following ISO C99 standard, Annex G,
+# more precisely Section 7.26 "Future library directions"
+
+# Rule [conj]: log10(conj(z)) = conj(log10(z))
+
+# log10(nan + i*inf) = +inf + i*nan
+0 0 2 +inf 2 nan                   2 nan 2 +inf  N N
+0 0 2 +inf 2 nan                   2 nan 2 -inf  N N
+
+# log10(nan + i*nan) = nan + i*nan
+0 0 2 nan 2 nan                    2 nan 2 nan   N N
+
+# log10(nan + i*y) = nan + i*nan for finite y
+0 0 2 nan 2 nan                    2 nan 2  1    N N
+0 0 2 nan 2 nan                    2 nan 2 +0    N N
+0 0 2 nan 2 nan                    2 nan 2 -0    N N
+0 0 2 nan 2 nan                    2 nan 2 -1    N N
+
+# log10(+/-inf + i*nan) = +inf + i*nan
+0 0 2 +inf 2 nan                   2 +inf 2 nan  N N
+0 0 2 +inf 2 nan                   2 -inf 2 nan  N N
+
+# log10(x + i*nan) = nan+i*nan for finite x
+0 0 2 nan 2 nan                    2  1 2 nan    N N
+0 0 2 nan 2 nan                    2 +0 2 nan    N N
+0 0 2 nan 2 nan                    2 -0 2 nan    N N
+0 0 2 nan 2 nan                    2 -1 2 nan    N N
+
+# log10(x + i*inf) = +inf + i*pi/2/log(10) for finite x
+0 + 2 +inf 53 +0x15d47c4cb2fba1p-53  2  1 2 +inf   N N
+0 + 2 +inf 53 +0x15d47c4cb2fba1p-53  2 +0 2 +inf   N N
+0 + 2 +inf 53 +0x15d47c4cb2fba1p-53  2 -0 2 +inf   N N
+0 + 2 +inf 53 +0x15d47c4cb2fba1p-53  2 -1 2 +inf   N N
+# by [conj]: log10(x - i*inf) = +inf - i*pi/2/log(10) for finite x
+0 - 2 +inf 53 -0x15d47c4cb2fba1p-53  2  1 2 -inf  N N
+0 - 2 +inf 53 -0x15d47c4cb2fba1p-53  2 +0 2 -inf  N N
+0 - 2 +inf 53 -0x15d47c4cb2fba1p-53  2 -0 2 -inf  N N
+0 - 2 +inf 53 -0x15d47c4cb2fba1p-53  2 -1 2 -inf  N N
+
+# log10(-inf + i*y) = +inf + i*pi/log(10) for finite positive-signed y
+0 + 2 +inf 53 0x15d47c4cb2fba1p-52  2 -inf 2  1   N N
+0 + 2 +inf 53 0x15d47c4cb2fba1p-52  2 -inf 2 +0   N N
+# by [conj]: log10(-inf+i*y) = +inf - i*pi/log(10) for finite negative-signed y
+0 - 2 +inf 53 -0x15d47c4cb2fba1p-52  2 -inf 2 -1  N N
+0 - 2 +inf 53 -0x15d47c4cb2fba1p-52  2 -inf 2 -0  N N
+
+# log10(+inf + i*y) = +inf + i*0 for finite positive-signed y
+0 0 2 +inf 2 +0                    2 +inf 2  1   N N
+0 0 2 +inf 2 +0                    2 +inf 2 +0   N N
+# by [conj]: log10(+inf + i*y) = +inf - i*0 for finite negative-signed y
+0 0 2 +inf 2 -0                    2 +inf 2 -1   N N
+0 0 2 +inf 2 -0                    2 +inf 2 -0   N N
+
+# log10(-inf + i*inf) = +inf + i*(3*pi/4)/log(10)
+0 + 2 +inf 53 0x105f5d39863cb9p-52  2 -inf 2 +inf N N
+# by [conj]: log10(-inf - i*inf) = +inf - i*(3*pi/4)/log(10)
+0 - 2 +inf 53 -0x105f5d39863cb9p-52 2 -inf 2 -inf N N
+
+# log10(+inf + i*inf) = +inf + i*(pi/4)/log(10)
+0 + 2 +inf 53 0x15d47c4cb2fba1p-54  2 +inf 2 +inf N N
+# by [conj]: log10(+inf - i*inf) = +inf - i*(pi/4)/log(10)
+0 - 2 +inf 53 -0x15d47c4cb2fba1p-54 2 +inf 2 -inf N N
+
+# log10(-0 + i*0) = -inf + i*pi/log(10)
+0 + 2 -inf 53 +0x15d47c4cb2fba1p-52   2 -0 2 +0     N N
+0 - 2 -inf 53 -0x15d47c4cb2fba1p-52   2 -0 2 -0     N N
+0 + 2 -inf 53 +0x15d47c4cb2fba1p-52   2 -0 2 +0     N U
+0 + 2 -inf 53 -0xaea3e26597ddp-47     2 -0 2 -0     N U
+0 - 2 -inf 53 +0xaea3e26597ddp-47     2 -0 2 +0     N D
+0 - 2 -inf 53 -0x15d47c4cb2fba1p-52   2 -0 2 -0     N D
+
+# log10(+0 + i*0) = -inf + i*0
+0 0 2 -inf 2 +0                      2 +0 2 +0    N N
+0 0 2 -inf 2 -0                      2 +0 2 -0    N N
+
+# log10(+1 +- i*0) = +0 +- i*0
+0 0 2 +0 2 +0                      2 1 2 +0    N N
+0 0 2 +0 2 -0                      2 1 2 -0    N N
+
+# log10(10 +- i*0) = 1 +- i*0
+0 0 2 1 2 +0                      4 10 2 +0    N N
+0 0 2 1 2 -0                      4 10 2 -0    N N
+
+# log10(100 +- i*0) = 1 +- i*0
+0 0 2 2 2 +0                      5 100 2 +0    N N
+0 0 2 2 2 -0                      5 100 2 -0    N N
+
+# log10(-1 +- i*0) = +0 +- i*pi/log(10)
+0 + 2 +0 53 +0x15d47c4cb2fba1p-52   2 -1 2 +0     N N
+0 - 2 +0 53 -0x15d47c4cb2fba1p-52   2 -1 2 -0     N N
+
+# log10(x + i*y) with either x or y zero and the other non-zero
+- 0 53 0x13afeb354b7d97p-52  2                   0     5  0x11 2     0 N N
+- + 53 0x13afeb354b7d97p-52 53  0x15d47c4cb2fba1p-53   2     0 5  0x11 N N
+- + 53 0x1475c655fbc11p-48  53  0x15d47c4cb2fba1p-52   5 -0x13 2    +0 N N
+- - 53 0x1475c655fbc11p-48  53 -0x15d47c4cb2fba1p-52   5 -0x13 2    -0 N N
+- - 53 0x1475c655fbc11p-48  53 -0x15d47c4cb2fba1p-53   2     0 5 -0x13 N N
+
+- + 53 0x15c9a3209bf97fp-52 53  0x15d47c4cb2fba1p-52   5 -0x17 2    +0 Z U
+- - 53 0x15c9a3209bf97fp-52 53 -0x15d47c4cb2fba1p-52   5 -0x17 2    -0 N D
++ - 53 0x2b93464137f3p-45   53  0xaea3e26597ddp-47     5 -0x17 2    +0 U Z
+- - 53 0x15c9a3209bf97fp-52 53 -0x15d47c4cb2fba1p-52   5 -0x17 2    -0 D N
+- - 53 0x1d1cda1a0c996dp-52 53  0xaea3e26597ddp-47     7 -0x42 2    +0 Z D
++ + 53 0xe8e6d0d064cb7p-51  53 -0xaea3e26597ddp-47     7 -0x42 2    -0 N Z
++ + 53 0xe8e6d0d064cb7p-51  53  0x15d47c4cb2fba1p-52   7 -0x42 2    +0 U N
+- + 53 0x1d1cda1a0c996dp-52 53 -0xaea3e26597ddp-47     7 -0x42 2    -0 D U
+
+# huge values
++ 0 53 0x134413509f79ffp-44 2   0                   2 0x1p1024 2    +0 U U
++ 0 53 0x134413509f79ffp-43 2   0                   2 0x1p2048 2    +0 U U
++ 0 53 0x134413509f79ffp-42 2   0                   2 0x1p4096 2    +0 U U
++ 0 53 0x134413509f79ffp-41 2   0                   2 0x1p8192 2    +0 U U
++ 0 53 0x134413509f79ffp-40 2   0                   2 0x1p16384 2   +0 U U
++ 0 53 0x134413509f79ffp-39 2   0                   2 0x1p32768 2   +0 U U
++ 0 53 0x134413509f79ffp-38 2   0                   2 0x1p65536 2   +0 U U
++ 0 53 0x134413509f79ffp-37 2   0                   2 0x1p131072 2  +0 U U
++ 0 53 0x134413509f79ffp-36 2   0                   2 0x1p262144 2  +0 U U
++ 0 53 0x134413509f79ffp-35 2   0                   2 0x1p524288 2  +0 U U
++ 0 53 0x134413509f79ffp-34 2   0                   2 0x1p1048576 2 +0 U U
++ 0 53 0x134413509f79ffp-33 2   0                   2 0x1p2097152 2 +0 U U
++ 0 53 0x134413509f79ffp-32 2   0                   2 0x1p4194304 2 +0 U U
++ 0 53 0x134413509f79ffp-31 2   0                   2 0x1p8388608 2 +0 U U
++ 0 53 0x134413509f79ffp-30 2   0                  2 0x1p16777216 2 +0 U U
++ 0 53 0x134413509f79ffp-29 2   0                  2 0x1p33554432 2 +0 U U
++ 0 53 0x134413509f79ffp-28 2   0                  2 0x1p67108864 2 +0 U U
++ 0 53 0x134413509f79ffp-27 2   0                 2 0x1p134217728 2 +0 U U
++ 0 53 0x134413509f79ffp-26 2   0                 2 0x1p268435456 2 +0 U U
++ 0 53 0x134413509f79ffp-25 2   0                 2 0x1p536870912 2 +0 U U
++ + 53 0x13467bd3098defp-44 53 0x15d47c4cb2fba1p-54 2 0x1p1024 2 0x1p1024 U U
++ + 53 0x13454791d483f7p-43 53 0x15d47c4cb2fba1p-54 2 0x1p2048 2 0x1p2048 U U
++ + 53 0x1344ad7139fefbp-42 53 0x15d47c4cb2fba1p-54 2 0x1p4096 2 0x1p4096 U U
++ + 53 0x13446060ecbc7dp-41 53 0x15d47c4cb2fba1p-54 2 0x1p8192 2 0x1p8192 U U
++ + 53 0x9a21cec630d9fp-39  53 0x15d47c4cb2fba1p-54 2 0x1p16384 2 0x1p16384 U U
++ + 53 0x13442694b2ca9fp-39 53 0x15d47c4cb2fba1p-54 2 0x1p32768 2 0x1p32768 U U
++ + 53 0x13441cf2a9224fp-38 53 0x15d47c4cb2fba1p-54 2 0x1p65536 2 0x1p65536 U U
++ + 53 0x13441821a44e27p-37  2  0x3p-3            2 0x1p131072 2 0x1p131072 U U
++ + 53 0x134415b921e413p-36  2  0x3p-3            2 0x1p262144 2 0x1p262144 U U
++ + 53 0x13441484e0af09p-35  2  0x3p-3            2 0x1p524288 2 0x1p524288 U U
++ + 53 0x4d104fab00521p-32   2  0x3p-3          2 0x1p1048576 2 0x1p1048576 U U
++ + 53 0x9a209ced7e3a1p-32   2  0x3p-3          2 0x1p2097152 2 0x1p2097152 U U
++ + 53 0x1344137727a0a1p-32  2  0x3p-3          2 0x1p4194304 2 0x1p4194304 U U
++ + 53 0x13441363e38d5p-27   2  0x3p-3          2 0x1p8388608 2 0x1p8388608 U U
++ + 53 0x268826b483075p-27   2  0x3p-3        2 0x1p16777216 2 0x1p16777216 U U
++ + 53 0x4d104d55c1fb5p-27   2  0x3p-3        2 0x1p33554432 2 0x1p33554432 U U
++ + 53 0x9a209a983fe35p-27   2  0x3p-3        2 0x1p67108864 2 0x1p67108864 U U
++ + 53 0x4d104d474eecdp-25   2  0x3p-3      2 0x1p134217728 2 0x1p134217728 U U
++ + 53 0x9a209a89ccd4dp-25   2  0x3p-3      2 0x1p268435456 2 0x1p268435456 U U
+# Due to intermediate overflow, the following result has wrong real part
+# 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 0x13441350ec8a4dp-25  2  0x3p-3    2 0x1p536870912 2 0x1p536870912 U U
diff --git a/tests/tlog10.c b/tests/tlog10.c
new file mode 100644
index 0000000..5591b06
--- /dev/null
+++ b/tests/tlog10.c
@@ -0,0 +1,37 @@
+/* tlog10 -- test file for mpc_log10.
+
+Copyright (C) 2012 INRIA
+
+This file is part of GNU MPC.
+
+GNU MPC is free software; you can redistribute it and/or modify it under
+the terms of the GNU Lesser General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
+more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program. If not, see http://www.gnu.org/licenses/ .
+*/
+
+#include <stdlib.h>
+#include "mpc-tests.h"
+
+int
+main (void)
+{
+  DECL_FUNC (CC, f, mpc_log10);
+
+  test_start ();
+
+  data_check (f, "log10.dat");
+  tgeneric (f, 2, 512, 7, 128);
+
+  test_end ();
+
+  return 0;
+}