Synchronize with r1413 from trunk.benchs_tests

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

12 files changed, 134 insertions, 22 deletions

diff --git a/NEWS b/NEWS
index b754788..4ef1540 100644
--- a/NEWS
+++ b/NEWS
@@ -1,11 +1,13 @@
 Recent changes in the trunk:
   - Minimally required library version: mpfr 3.0.0
+  - Improved speed for corner cases of mpc_asin, mpc_sin, see
+    http://lists.gforge.inria.fr/pipermail/mpc-discuss/2013-December/001266.html
 
 Changes in version 1.0.2:
-   - Fixed mpc_atan, mpc_atanh for (+-0, +-1), see
-     http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57994#c7
-   - Fixed mpc_log10 for purely imaginary argument, see
-     http://lists.gforge.inria.fr/pipermail/mpc-discuss/2012-September/001208.html
+  - Fixed mpc_atan, mpc_atanh for (+-0, +-1), see
+    http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57994#c7
+  - Fixed mpc_log10 for purely imaginary argument, see
+    http://lists.gforge.inria.fr/pipermail/mpc-discuss/2012-September/001208.html
 
 Changes in version 1.0.1:
   - Switched to automake 1.11.6, see
@@ -149,7 +151,7 @@ Changes in version 0.5 ("Aconitum neomontanum"):
     - mpc_sqrt with directed rounding
 
 
-Copyright (C) 2008, 2009, 2010, 2011, 2012 INRIA
+Copyright (C) 2008, 2009, 2010, 2011, 2012, 2013, 2014 INRIA
 
 Copying and distribution of this file, with or without modification,
 are permitted in any medium without royalty provided the copyright
diff --git a/TODO b/TODO
index 4a7f47f..06d041c 100644
--- a/TODO
+++ b/TODO
@@ -1,3 +1,7 @@
+From Karim Belabas 9 Jan 2014:
+Implement Hurwitz(s,x) -> gives Zeta for x=1.
+Cf http://arxiv.org/abs/1309.2877
+
 From Andreas Enge 27 August 2012:
 Implement im(atan(x+i*y)) as
 1/4 * [log1p (4y / (x^2 +(1-y)^2))]
@@ -33,6 +37,7 @@ sage: %timeit asin(MPComplexField()(1,1))
 10000 loops, best of 3: 83.7 us per loop
 sage: %timeit asin(MPComplexField()(1,1e-1000))
 100 loops, best of 3: 17 ms per loop
+-> should be much faster with revision 1402 (check)
 
 Same for acos:
 sage: %timeit acos(MPComplexField()(1,1))      
diff --git a/doc/algorithms.bib b/doc/algorithms.bib
index 2232a0e..eb4855e 100644
--- a/doc/algorithms.bib
+++ b/doc/algorithms.bib
@@ -36,3 +36,11 @@
   number =	 4,
   pages =	 {359--367}
 }
+
+@Misc{CrTh10,
+  author = 	 {John E. Cremona and Thotsaphon Thongjunthug},
+  title = 	 {The complex {AGM}, periods of elliptic curves over {\C} and complex elliptic logarithms},
+  howpublished = {\url{http://arxiv.org/abs/1011.0914}},
+  year = 	 2010,
+  note = 	 {31 pages}}
+
diff --git a/doc/algorithms.tex b/doc/algorithms.tex
index 669d651..d9acd15 100644
--- a/doc/algorithms.tex
+++ b/doc/algorithms.tex
@@ -17,6 +17,7 @@
 \newcommand {\Ulp}{{\operatorname {ulp}}}
 \DeclareMathOperator{\Exp}{\operatorname {Exp}}
 \newcommand {\atantwo}{\operatorname {atan2}}
+\newcommand {\asin}{\operatorname {asin}}
 \newcommand{\error}{\operatorname {error}}
 \newcommand{\relerror}{\operatorname {relerror}}
 \newcommand{\Norm}{\operatorname {N}}
@@ -1142,6 +1143,11 @@ k_R=\left\{
 \right.
 \end{equation*}
 
+\subsection {\texttt {mpc\_asin}}
+
+For $0 \leq y \leq 1$, we have $|\Re \asin (1 \pm iy) - \pi/2|
+\leq y^{1/2}$, and $|\Im \asin (1 \pm iy) \mp y^{1/2}| \leq (1/12) y^{3/2}$.
+
 \subsection {\texttt {mpc\_pow}}
 
 The main issue for the power function is to be able to recognize when the
diff --git a/m4/mpc.m4 b/m4/mpc.m4
index 414229b..1204c3b 100644
--- a/m4/mpc.m4
+++ b/m4/mpc.m4
@@ -90,7 +90,7 @@ AC_DEFUN([MPC_C_CHECK_FLAG], [
 AC_DEFUN([MPC_C_CHECK_WARNINGCFLAGS], [
   AC_REQUIRE([AC_PROG_GREP])
   if echo $VERSION | grep -c dev >/dev/null 2>&1 ; then
-    if test "x$GCC" = "xyes" -a "x$compiler" != "xicc" -a "x$compiler" != "xg++"; then
+    if test "x$GCC" = "xyes" -a "x$compiler" != "xicc"; then
       # enable -Werror for myself (Andreas Enge)
       if test "x$USER" = "xenge"; then
          MPC_C_CHECK_FLAG(-Werror)
@@ -107,6 +107,7 @@ AC_DEFUN([MPC_C_CHECK_WARNINGCFLAGS], [
       MPC_C_CHECK_FLAG(-Wstrict-prototypes)
       MPC_C_CHECK_FLAG(-Wmissing-prototypes)
       MPC_C_CHECK_FLAG(-Wno-unused-value)
+      MPC_C_CHECK_FLAG(-Wlogical-op)
     fi
   fi
 ])
diff --git a/src/acos.c b/src/acos.c
index 6117786..b95387d 100644
--- a/src/acos.c
+++ b/src/acos.c
@@ -24,7 +24,7 @@ along with this program. If not, see http://www.gnu.org/licenses/ .
 int
 mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
 {
-  int inex_re, inex_im, inex;
+  int inex_re, inex_im, inex, loop = 0;
   mpfr_prec_t p_re, p_im, p;
   mpc_t z1;
   mpfr_t pi_over_2;
@@ -190,7 +190,8 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
   mpfr_init2 (pi_over_2, p);
   for (;;)
     {
-      p += mpc_ceil_log2 (p) + 3;
+      loop ++;
+      p += (loop <= 2) ? mpc_ceil_log2 (p) + 3 : p / 2;
 
       mpfr_set_prec (mpc_realref(z1), p);
       mpfr_set_prec (pi_over_2, p);
diff --git a/src/asin.c b/src/asin.c
index 2d18489..a92421a 100644
--- a/src/asin.c
+++ b/src/asin.c
@@ -1,6 +1,6 @@
 /* mpc_asin -- arcsine of a complex number.
 
-Copyright (C) 2009, 2010, 2011, 2012 INRIA
+Copyright (C) 2009, 2010, 2011, 2012, 2013 INRIA
 
 This file is part of GNU MPC.
 
@@ -20,13 +20,74 @@ along with this program. If not, see http://www.gnu.org/licenses/ .
 
 #include "mpc-impl.h"
 
+/* Special case op = 1 + i*y for tiny y (see algorithms.tex).
+   Return 0 if special formula fails, otherwise put in z1 the approximate
+   value which needs to be converted to rop.
+   z1 is a temporary variable with enough precision.
+ */
+static int
+mpc_asin_special (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd, mpc_ptr z1)
+{
+  mpfr_exp_t ey = mpfr_get_exp (mpc_imagref (op));
+  mpfr_t abs_y;
+  mpfr_prec_t p;
+  int inex;
+
+  /* |Re(asin(1+i*y)) - pi/2| <= y^(1/2) */
+  if (ey >= 0 || ((-ey) / 2 < mpfr_get_prec (mpc_realref (z1))))
+    return 0;
+  if (mpfr_get_prec (mpc_imagref (z1)) > (-ey) + 2) /* p > -ey + 2 */
+    return 0;
+
+  mpfr_const_pi (mpc_realref (z1), MPFR_RNDN);
+  mpfr_div_2exp (mpc_realref (z1), mpc_realref (z1), 1, MPFR_RNDN); /* exact */
+  p = mpfr_get_prec (mpc_realref (z1));
+  /* if z1 has precision p, the error on z1 is 1/2*ulp(z1) = 2^(-p) so far,
+     and since ey <= -2p, then y^(1/2) <= 1/2*ulp(z1) too, thus the total
+     error is bounded by ulp(z1) */
+  if (!mpfr_can_round (mpc_realref(z1), p, MPFR_RNDN, MPFR_RNDZ,
+                       mpfr_get_prec (mpc_realref(rop)) +
+                       (MPC_RND_RE(rnd) == MPFR_RNDN)))
+    return 0;
+  
+  /* |Im(asin(1+i*y)) - y^(1/2)| <= (1/12) * y^(3/2) for y >= 0 (err >= 0)
+     |Im(asin(1-i*y)) + y^(1/2)| <= (1/12) * y^(3/2) for y >= 0 (err <= 0) */
+  abs_y[0] = mpc_imagref (op)[0];
+  if (mpfr_signbit (mpc_imagref (op)))
+    MPFR_CHANGE_SIGN (abs_y);
+  inex = mpfr_sqrt (mpc_imagref (z1), abs_y, MPFR_RNDN); /* error <= 1/2 ulp */
+  if (mpfr_signbit (mpc_imagref (op)))
+    MPFR_CHANGE_SIGN (mpc_imagref (z1));
+  /* If z1 has precision p, the error on z1 is 1/2*ulp(z1) = 2^(-p) so far,
+     and (1/12) * y^(3/2) <= (1/8) * y * y^(1/2) <= 2^(ey-3)*2^p*ulp(y^(1/2))
+     thus for p+ey-3 <= -1 we have (1/12) * y^(3/2) <= (1/2) * ulp(y^(1/2)),
+     and the total error is bounded by ulp(z1).
+     Note: if y^(1/2) is exactly representable, and ends with many zeroes,
+     then mpfr_can_round below will fail; however in that case we know that
+     Im(asin(1+i*y)) is away from +/-y^(1/2) wrt zero. */
+  if (inex == 0) /* enlarge im(z1) so that the inexact flag is correct */
+    {
+      if (mpfr_signbit (mpc_imagref (op)))
+        mpfr_nextbelow (mpc_imagref (z1));
+      else
+        mpfr_nextabove (mpc_imagref (z1));
+      return 1;
+    }
+  p = mpfr_get_prec (mpc_imagref (z1));
+  if (!mpfr_can_round (mpc_imagref(z1), p - 1, MPFR_RNDA, MPFR_RNDZ,
+                      mpfr_get_prec (mpc_imagref(rop)) +
+                      (MPC_RND_IM(rnd) == MPFR_RNDN)))
+    return 0;
+  return 1;
+}
+
 int
 mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
 {
-  mpfr_prec_t p, p_re, p_im, incr_p = 0;
+  mpfr_prec_t p, p_re, p_im;
   mpfr_rnd_t rnd_re, rnd_im;
   mpc_t z1;
-  int inex;
+  int inex, loop = 0;
 
   /* special values */
   if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op)))
@@ -148,11 +209,15 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
   {
     mpfr_exp_t ex, ey, err;
 
-    p += mpc_ceil_log2 (p) + 3 + incr_p; /* incr_p is zero initially */
-    incr_p = p / 2;
+    loop ++;
+    p += (loop <= 2) ? mpc_ceil_log2 (p) + 3 : p / 2;
     mpfr_set_prec (mpc_realref(z1), p);
     mpfr_set_prec (mpc_imagref(z1), p);
 
+    /* try special code for 1+i*y with tiny y */
+    if (loop == 1 && mpc_asin_special (rop, op, rnd, z1))
+      break;
+
     /* z1 <- z^2 */
     mpc_sqr (z1, op, MPC_RNDNN);
     /* err(x) <= 1/2 ulp(x), err(y) <= 1/2 ulp(y) */
diff --git a/src/sin_cos.c b/src/sin_cos.c
index 0bb00cf..7051708 100644
--- a/src/sin_cos.c
+++ b/src/sin_cos.c
@@ -299,11 +299,28 @@ mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
       mpfr_t s, c, sh, ch, sch, csh;
       mpfr_prec_t prec;
       int ok;
-      int inex_re, inex_im, inex_sin, inex_cos;
+      int inex_re, inex_im, inex_sin, inex_cos, loop = 0;
 
       prec = 2;
       if (rop_sin != NULL)
-         prec = MPC_MAX (prec, MPC_MAX_PREC (rop_sin));
+        {
+          mp_prec_t er, ei;
+          prec = MPC_MAX (prec, MPC_MAX_PREC (rop_sin));
+          /* since the Taylor expansion of sin(x) at x=0 is x - x^3/6 + O(x^5),
+             if x <= 2^(-p), then the second term/x is about 2^(-2p)/6, thus we
+             need at least 2p+3 bits of precision. This is true only when x is
+             exactly representable in the target precision. */
+          if (MPC_MAX_PREC (op) <= prec)
+            {
+              er = mpfr_get_exp (mpc_realref (op));
+              ei = mpfr_get_exp (mpc_imagref (op));
+              /* consider the maximal exponent only */
+              er = (er < ei) ? ei : er;
+              if (er < 0)
+                if (prec < 2 * (mp_prec_t) (-er) + 3)
+                  prec = 2 * (mp_prec_t) (-er) + 3;
+            }
+        }
       if (rop_cos != NULL)
          prec = MPC_MAX (prec, MPC_MAX_PREC (rop_cos));
 
@@ -315,8 +332,9 @@ mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
       mpfr_init2 (csh, 2);
 
       do {
+         loop ++;
          ok = 1;
-         prec += mpc_ceil_log2 (prec) + 5;
+         prec += (loop <= 2) ? mpc_ceil_log2 (prec) + 5 : prec / 2;
 
          mpfr_set_prec (s, prec);
          mpfr_set_prec (c, prec);
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 77e5866..61e5583 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -25,8 +25,9 @@ AM_LDFLAGS = -no-install
 # LOADLIBES (documented in the "GNU make" manual and equivalent to LDLIBS)
 # enables to compile a program foo.c in the test directory by simply doing
 # "make foo".
-LOADLIBES=$(DEFS) -I$(top_srcdir)/src -I$(top_builddir) $(CPPFLAGS)	\
-  $(CFLAGS) -L$(top_builddir)/tests/.libs -lmpc-tests -lmpc -lm $(LIBS)
+LOADLIBES=$(DEFS) $(AM_CPPFLAGS) $(CPPFLAGS) $(CFLAGS)	\
+  $(top_builddir)/tests/.libs/libmpc-tests.a		\
+  $(top_builddir)/src/.libs/libmpc.a $(LIBS)
 
 check_PROGRAMS = tabs tacos tacosh tadd tadd_fr tadd_si tadd_ui targ	\
   tasin tasinh tatan tatanh tconj tcos tcosh tdiv tdiv_2si tdiv_2ui	\
diff --git a/tests/asin.dat b/tests/asin.dat
index 9793fcc..e93e8c5 100644
--- a/tests/asin.dat
+++ b/tests/asin.dat
@@ -124,3 +124,5 @@
 + - 53 0x189BF9EC7FCD5Bp-54 53 0x1206ECFA94614Bp-50 53 17 53 42 N N
 - + 2 1.5 2 6 2 96 2 0x1p-8 N N
 - - 8 0xC9p-7 8 0x15p-2 2 96 2 0x1p-8 N N
+- - 53 0x3243f6a8885a3p-49 53 0x1p-500 53 1 53 0x1p-1000 N N
+- + 53 0x3243f6a8885a3p-49 53 -0x1p-500 53 1 53 -0x1p-1000 N N
diff --git a/tests/tio_str.c b/tests/tio_str.c
index b7daf9d..b4e0c34 100644
--- a/tests/tio_str.c
+++ b/tests/tio_str.c
@@ -1,6 +1,6 @@
 /* tio_str-- Test file for mpc_inp_str and mpc_out_str.
 
-Copyright (C) 2009, 2011, 2013 INRIA
+Copyright (C) 2009, 2011, 2013, 2014 INRIA
 
 This file is part of GNU MPC.
 
@@ -21,9 +21,9 @@ along with this program. If not, see http://www.gnu.org/licenses/ .
 #include "mpc-tests.h"
 
 #ifdef HAVE_UNISTD_H
-#if !defined _POSIX_C_SOURCE
+#ifndef _POSIX_C_SOURCE
 #define _POSIX_C_SOURCE 1 /* apparently needed on Darwin */
-#endif /* _POSIX_C_SOURCE */
+#endif
 #include <unistd.h> /* for dup, dup2, STDIN_FILENO and STDOUT_FILENO */
 #else
 #define STDIN_FILENO 0
diff --git a/tests/tpow_ld.c b/tests/tpow_ld.c
index a02ea6e..a4e3ab9 100644
--- a/tests/tpow_ld.c
+++ b/tests/tpow_ld.c
@@ -1,6 +1,6 @@
 /* tpow_ld -- test file for mpc_pow_ld.
 
-Copyright (C) 2009 INRIA
+Copyright (C) 2009, 2014 INRIA
 
 This file is part of GNU MPC.
 
@@ -33,6 +33,9 @@ main (void)
   if (mpc_cmp_si_si (z, -9, 46) != 0)
     {
       printf ("Error for mpc_pow_ld (1)\n");
+      printf ("expected (-9 46)\ngot ");
+      mpc_out_str (stdout, 10, 0, z, MPC_RNDNN);
+      printf ("\n");
       exit (1);
     }
   mpc_clear (z);