Merged the latest changes (in particular the fix in sub1.c) from the

trunk.

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10388 280ebfd0-de03-0410-8827-d642c229c3f4

14 files changed, 357 insertions, 91 deletions

diff --git a/TODO b/TODO
index cda9815fb..700cc5875 100644
--- a/TODO
+++ b/TODO
@@ -269,6 +269,9 @@ Table of contents:
 5. Efficiency
 ##############################################################################
 
+- Fredrik Johansson reports that mpfr_ai is slow for large arguments: an
+  asymptotic expansion should be used (once done, remove REDUCE_EMAX from
+  tests/tai.c and update the description in mpfr.texi).
 - for exp(x), Fredrik Johansson reports a 20% speed improvement starting from
   4000 bits, and up to a 75% memory improvement in his Arb implementation, by
   using recursive instead of iterative binary splitting:
diff --git a/acinclude.m4 b/acinclude.m4
index 76e3e5a79..be917becb 100644
--- a/acinclude.m4
+++ b/acinclude.m4
@@ -209,7 +209,11 @@ fi
 dnl Check for attribute constructor and destructor
 MPFR_CHECK_CONSTRUCTOR_ATTR()
 
-dnl Check for POSIX Thread.
+dnl Check for POSIX Thread. Since the AX_PTHREAD macro is not standard
+dnl (it is provided by autoconf-archive), we need to detect whether it
+dnl is left unexpanded, otherwise the configure script won't fail and
+dnl "make distcheck" won't give any error, yielding buggy tarballs!
+m4_pattern_forbid([AX_PTHREAD])
 AX_PTHREAD([])
 
 dnl Check for ISO C11 Thread
diff --git a/doc/README.dev b/doc/README.dev
index 458017eea..4d9780df0 100644
--- a/doc/README.dev
+++ b/doc/README.dev
@@ -27,9 +27,9 @@ need some GNU development utilities: aclocal, autoheader, automake,
 autoconf 2.60 (at least), libtoolize, and the AX_PTHREAD macro from
 <http://www.gnu.org/software/autoconf-archive/ax_pthread.html>.
 As some files like "configure" are not part of the Subversion
-repository, you first need to run "autoreconf -i" (or autogen.sh,
+repository, you first need to run "autoreconf -i" (or ./autogen.sh,
 which could be used later to update the config files). Then you can
-run "configure" in the usual way (see the INSTALL file, but note that
+run ./configure in the usual way (see the INSTALL file, but note that
 there are no patches to apply, and the URLs are not valid since the
 corresponding version has not been released yet).
 
diff --git a/src/log_ui.c b/src/log_ui.c
index 3722dad0e..b544e0d64 100644
--- a/src/log_ui.c
+++ b/src/log_ui.c
@@ -24,53 +24,63 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
 #include "mpfr-impl.h"
 
 /* FIXME: mpfr_log_ui is much slower than mpfr_log on some values of n,
-   e.g. something like 8 times as slow for n around ULONG_MAX/3 on an
-   x86_64 Linux machine. */
-
-/* Auxiliary function: Compute using binary splitting the sum
-   -sum((-x)^(n+1-n1)/n, n = n1..n2-1) where x = p/2^k,
-   -1/3 <= x <= 1/3. For n1=1 we get -sum((-x)^n/n, n = 1..n2-1).
-   Numerator is P[0], denominator is Q[0],
-   Need 1+ceil(log(n2-n1)/log(2)) cells in P[],Q[].
+   e.g. about 4 times as slow for n around ULONG_MAX/3 on an
+   x86_64 Linux machine, for 10^6 bits of precision. The reason is that
+   for say n=6148914691236517205 and prec=10^6, the value of T computed
+   has more than 50M bits, which is much more than needed. Indeed the
+   binary splitting algorithm for series with a finite radius of convergence
+   gives rationals of size n*log(n) for a target precision n. One might
+   truncate the rationals inside the algorithm, but then the error analysis
+   should be redone. */
+
+/* Cf http://www.ginac.de/CLN/binsplit.pdf: the Taylor series of log(1+x)
+   up to order N for x=p/2^k is T/(B*Q).
+   P[0] <- (-p)^(n2-n1) [with opposite sign when n1=1]
+   q <- k*(n2-n1) [corresponding to Q[0] = 2^q]
+   B[0] <- n1 * (n1+1) * ... * (n2-1)
+   T[0] <- B[0]*Q[0] * S(n1,n2)
+   where S(n1,n2) = -sum((-x)^(i-n1+1)/i, i=n1..n2-1)
+   Assumes p is odd or zero, and -1/3 <= x = p/2^k <= 1/3.
 */
 static void
-S (mpz_t *P, mpz_t *Q, unsigned long n1, unsigned long n2,
-   long p, unsigned long k)
+S (mpz_t *P, unsigned long *q, mpz_t *B, mpz_t *T, unsigned long n1,
+   unsigned long n2, long p, unsigned long k, int need_P)
 {
   MPFR_ASSERTD (n1 < n2);
+  MPFR_ASSERTD (p == 0 || ((unsigned long) p & 1) != 0);
   if (n2 == n1 + 1)
     {
-      mpz_set_si (P[0], p);
-      mpz_set_ui (Q[0], n1);
-      mpz_mul_2exp (Q[0], Q[0], k);
+      mpz_set_si (P[0], (n1 == 1) ? p : -p);
+      *q = k;
+      mpz_set_ui (B[0], n1);
+      /* T = B*Q*S where S = P/(B*Q) thus T = P */
+      mpz_set (T[0], P[0]);
+      /* since p is odd (or zero), there is no common factor 2 between
+         P and Q, or T and B */
     }
   else
     {
-      unsigned long m = (n1 / 2) + (n2 / 2) + (n1 & 1UL & n2);
+      unsigned long m = (n1 / 2) + (n2 / 2) + (n1 & 1UL & n2), q1;
       /* m = floor((n1+n2)/2) */
 
       MPFR_ASSERTD (n1 < m && m < n2);
-      S (P, Q, n1, m, p, k);
-      S (P + 1, Q + 1, m, n2, p, k);
-      /* P1/Q1 + (-p)^(m-n1)/2^(k(m-n1)) P2/Q2 =
-         (P1*2^(k(m-n1))*Q2 + Q1*(-p)^(m-n1)*P2)/(Q1*Q2*2^(k(m-n1)))
-         We know 2^(k(m-n1)) divides Q1, thus:
-         (P1*Q2 + (Q1/2^(k(m-n1)))*(-p)^(m-n1)*P2)/(Q1*Q2) */
-      mpz_mul (P[0], P[0], Q[1]);       /* P1*Q2 */
-      mpz_mul (Q[1], Q[0], Q[1]);       /* Q1*Q2 */
-      mpz_tdiv_q_2exp (Q[0], Q[0], k * (m - n1)); /* Q1/2^(k(m-n1)) */
-      mpz_mul (P[1], P[1], Q[0]);       /* Q1*P2/2^(k(m-n1)) */
-      mpz_swap (Q[0], Q[1]);            /* Q1*Q2 */
-      if (p > 0)
-        {
-          mpz_ui_pow_ui (Q[1], p, m - n1);       /* p^m */
-          if ((m - n1) & 1)
-            mpz_neg (Q[1], Q[1]);
-        }
-      else
-        mpz_ui_pow_ui (Q[1], -p, m - n1);       /* (-p)^m */
-      mpz_mul (P[1], P[1], Q[1]);       /* Q1*Q2*2^(km) */
-      mpz_add (P[0], P[0], P[1]);
+      S (P, q, B, T, n1, m, p, k, 1);
+      S (P + 1, &q1, B + 1, T + 1, m, n2, p, k, need_P);
+
+      /* T0 <- T0*B1*Q1 + P0*B0*T1 */
+      mpz_mul (T[1], T[1], P[0]);
+      mpz_mul (T[1], T[1], B[0]);
+      mpz_mul (T[0], T[0], B[1]);
+      /* Q[1] = 2^q1 */
+      mpz_mul_2exp (T[0], T[0], q1); /* mpz_mul (T[0], T[0], Q[1]) */
+      mpz_add (T[0], T[0], T[1]);
+      if (need_P)
+        mpz_mul (P[0], P[0], P[1]);
+      *q += q1; /* mpz_mul (Q[0], Q[0], Q[1]) */
+      mpz_mul (B[0], B[0], B[1]);
+
+      /* there should be no common factors 2 between P, Q and T,
+         since P is odd (or zero) */
     }
 }
 
@@ -79,10 +89,10 @@ mpfr_log_ui (mpfr_ptr x, unsigned long n, mpfr_rnd_t rnd_mode)
 {
   unsigned long k;
   mpfr_prec_t w; /* working precision */
-  mpz_t three_n, *P, *Q;
+  mpz_t three_n, *P, *B, *T;
   mpfr_t t, q;
   int inexact;
-  unsigned long N, lgN, i;
+  unsigned long N, lgN, i, kk;
   long p;
   MPFR_GROUP_DECL(group);
   MPFR_TMP_DECL(marker);
@@ -128,37 +138,58 @@ mpfr_log_ui (mpfr_ptr x, unsigned long n, mpfr_rnd_t rnd_mode)
   p = n > LONG_MAX ? - (long) - n : (long) n;
 
   MPFR_TMP_MARK(marker);
-  w = MPFR_PREC(x) + 10;
+  w = MPFR_PREC(x) + MPFR_INT_CEIL_LOG2 (MPFR_PREC(x)) + 10;
   MPFR_GROUP_INIT_2(group, w, t, q);
   MPFR_SAVE_EXPO_MARK (expo);
 
+  kk = k;
+  if (p != 0)
+    while ((p % 2) == 0) /* replace p/2^kk by (p/2)/2^(kk-1) */
+      {
+        p /= 2;
+        kk --;
+      }
+
   MPFR_ZIV_INIT (loop, w);
   for (;;)
     {
       mpfr_t tmp;
       unsigned int err;
-      /* we need at most w*log(2)/log(3) terms for an accuracy of w bits */
+      unsigned long q0;
+
+      /* we need at most w/log2(2^kk/|p|) terms for an accuracy of w bits */
       mpfr_init2 (tmp, 32);
-      /* 1354911329/2^31 is a 32-bit upper bound for log(2)/log(3) */
-      mpfr_set_ui_2exp (tmp, 1354911329, -31, MPFR_RNDU);
+      mpfr_set_ui (tmp, (p > 0) ? p : -p, MPFR_RNDU);
+      mpfr_log2 (tmp, tmp, MPFR_RNDU);
+      mpfr_ui_sub (tmp, kk, tmp, MPFR_RNDD);
       MPFR_ASSERTN (w <= ULONG_MAX);
-      mpfr_mul_ui (tmp, tmp, w, MPFR_RNDU);
+      mpfr_ui_div (tmp, w, tmp, MPFR_RNDU);
       N = mpfr_get_ui (tmp, MPFR_RNDU);
+      if (N < 2)
+        N = 2;
       lgN = MPFR_INT_CEIL_LOG2 (N) + 1;
       mpfr_clear (tmp);
-      P = (mpz_t *) MPFR_TMP_ALLOC (2 * lgN * sizeof (mpz_t));
-      Q = P + lgN;
+      P = (mpz_t *) MPFR_TMP_ALLOC (3 * lgN * sizeof (mpz_t));
+      B = P + lgN;
+      T = B + lgN;
       for (i = 0; i < lgN; i++)
         {
           mpz_init (P[i]);
-          mpz_init (Q[i]);
+          mpz_init (B[i]);
+          mpz_init (T[i]);
         }
-      S (P, Q, 1, N, p, k);
-      mpfr_set_z (t, P[0], MPFR_RNDN); /* t = P[0] * (1 + theta_1) */
-      mpfr_set_z (q, Q[0], MPFR_RNDN); /* q = Q[0] * (1 + theta_2) */
-      mpfr_div (t, t, q, MPFR_RNDN);   /* t = P[0]/Q[0] * (1 + theta_3)^3
+
+      S (P, &q0, B, T, 1, N, p, kk, 0);
+      /* mpz_mul (Q[0], B[0], Q[0]); */
+      /* mpz_mul_2exp (B[0], B[0], q0); */
+
+      mpfr_set_z (t, T[0], MPFR_RNDN); /* t = P[0] * (1 + theta_1) */
+      mpfr_set_z (q, B[0], MPFR_RNDN); /* q = B[0] * (1 + theta_2) */
+      mpfr_mul_2exp (q, q, q0, MPFR_RNDN); /* B[0]*Q[0] */
+      mpfr_div (t, t, q, MPFR_RNDN);   /* t = T[0]/(B[0]*Q[0])*(1 + theta_3)^3
                                             = log(n/2^k) * (1 + theta_4)^4
                                             for |theta_i| < 2^(-w) */
+
       /* argument reconstruction: add k*log(2) */
       mpfr_const_log2 (q, MPFR_RNDN);
       mpfr_mul_ui (q, q, k, MPFR_RNDN);
@@ -166,12 +197,17 @@ mpfr_log_ui (mpfr_ptr x, unsigned long n, mpfr_rnd_t rnd_mode)
       for (i = 0; i < lgN; i++)
         {
           mpz_clear (P[i]);
-          mpz_clear (Q[i]);
+          mpz_clear (B[i]);
+          mpz_clear (T[i]);
         }
-      /* the maximal error is 5 ulps for P/Q, since |(1+/-u)^4 - 1| < 5*u
+      /* The maximal error is 5 ulps for P/Q, since |(1+/-u)^4 - 1| < 5*u
          for u < 2^(-12), k ulps for k*log(2), and 1 ulp for the addition,
-         thus at most k+6 ulps */
-      err = MPFR_INT_CEIL_LOG2 (k + 6);
+         thus at most k+6 ulps.
+         Note that there might be some cancellation in the addition: the worst
+         case is when log(1 + p/2^kk) = log(2/3) ~ -0.405, and with n=3 which
+         gives k=2, thus we add 2*log(2) = 1.386. Thus in the worst case we
+         have an exponent decrease of 1, which accounts for +1 in the error. */
+      err = MPFR_INT_CEIL_LOG2 (k + 6) + 1;
       if (MPFR_LIKELY (MPFR_CAN_ROUND (t, w - err, MPFR_PREC(x), rnd_mode)))
         break;
 
diff --git a/src/mpfr-gmp.h b/src/mpfr-gmp.h
index 4f4f6901d..df437cff1 100644
--- a/src/mpfr-gmp.h
+++ b/src/mpfr-gmp.h
@@ -442,6 +442,21 @@ typedef const mp_limb_t *mpfr_limb_srcptr;
    float endianness, this is true everywhere we know of and it'd be a fairly
    strange system that did anything else.  */
 
+/* Note for MPFR: building with "gcc -std=c89 -pedantic -pedantic-errors"
+   fails if the bit-field type is unsigned long:
+
+     error: type of bit-field '...' is a GCC extension [-Wpedantic]
+
+   Though with -std=c99 no errors are obtained, this is still an extension
+   in C99, which says:
+
+     A bit-field shall have a type that is a qualified or unqualified version
+     of _Bool, signed int, unsigned int, or some other implementation-defined
+     type.
+
+   So, unsigned int should be better.
+*/
+
 #ifndef _MPFR_IEEE_FLOATS
 
 #ifdef HAVE_DOUBLE_IEEE_LITTLE_ENDIAN
@@ -450,10 +465,10 @@ union mpfr_ieee_double_extract
 {
   struct
     {
-      unsigned long manl:32;
-      unsigned long manh:20;
-      unsigned long exp:11;
-      unsigned long sig:1;
+      unsigned int manl:32;
+      unsigned int manh:20;
+      unsigned int exp:11;
+      unsigned int sig:1;
     } s;
   double d;
 };
@@ -465,10 +480,10 @@ union mpfr_ieee_double_extract
 {
   struct
     {
-      unsigned long sig:1;
-      unsigned long exp:11;
-      unsigned long manh:20;
-      unsigned long manl:32;
+      unsigned int sig:1;
+      unsigned int exp:11;
+      unsigned int manh:20;
+      unsigned int manl:32;
     } s;
   double d;
 };
diff --git a/src/mpfr-impl.h b/src/mpfr-impl.h
index ff1e351b5..e0167eff6 100644
--- a/src/mpfr-impl.h
+++ b/src/mpfr-impl.h
@@ -1148,7 +1148,7 @@ typedef union { mp_size_t s; mp_limb_t l; } mpfr_size_limb_t;
   MPFR_CACHE_ATTR mpfr_cache_t _cache = {{                           \
       {{ 0, MPFR_SIGN_POS, 0, (mp_limb_t *) 0 }}, 0, _func           \
       MPFR_DEFERRED_INIT_SLAVE_VALUE(_func)                          \
-    }};
+    }}
 
 /******************************************************
  ***************  Threshold parameters  ***************
diff --git a/src/sub1.c b/src/sub1.c
index d9ced1646..8733fe4dc 100644
--- a/src/sub1.c
+++ b/src/sub1.c
@@ -120,16 +120,15 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
          but we need to add a test first for complete coverage. */
       MPFR_SET_EXP (a, exp_b);
       MPFR_RNDRAW_EVEN (inexact, a, MPFR_MANT (b), bq,
-                        rnd_mode, MPFR_SIGN (a),
-                        if (MPFR_UNLIKELY (++ MPFR_EXP (a) > __gmpfr_emax))
-                        inexact = mpfr_overflow (a, rnd_mode, MPFR_SIGN (a)));
-      /* inexact = mpfr_set4 (a, b, rnd_mode, MPFR_SIGN (a));  */
+                        rnd_mode, MPFR_SIGN (a), ++ MPFR_EXP (a));
       if (inexact == 0)
         {
           /* a = b (Exact)
              But we know it isn't (Since we have to remove `c')
              So if we round to Zero, we have to remove one ulp.
              Otherwise the result is correctly rounded. */
+          /* An overflow is not possible. */
+          MPFR_ASSERTD (MPFR_EXP (a) <= __gmpfr_emax);
           if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (a)))
             {
               mpfr_nexttozero (a);
@@ -160,9 +159,14 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
              i.e. inexact= MPFR_EVEN_INEX */
           if (MPFR_UNLIKELY (inexact == MPFR_EVEN_INEX * MPFR_INT_SIGN (a)))
             {
-              mpfr_nexttozero (a);
+              if (MPFR_UNLIKELY (MPFR_EXP (a) > __gmpfr_emax))
+                mpfr_setmax (a, __gmpfr_emax);
+              else
+                mpfr_nexttozero (a);
               inexact = -MPFR_INT_SIGN (a);
             }
+          else if (MPFR_UNLIKELY (MPFR_EXP (a) > __gmpfr_emax))
+            inexact = mpfr_overflow (a, rnd_mode, MPFR_SIGN (a));
           MPFR_RET (inexact);
         }
     }
diff --git a/tests/memory.c b/tests/memory.c
index 64038e522..ac478117d 100644
--- a/tests/memory.c
+++ b/tests/memory.c
@@ -103,6 +103,9 @@ tests_memory_valid (void *ptr)
 }
 */
 
+/* FIXME: we might have an environment variable MPFR_TESTS_NO_MEMORY_LIMIT
+   to disable the memory limit check, for example if we want to compute
+   zillions of digits of Pi with ./tconst_pi 1000000000000 */
 static void
 tests_addsize (size_t size)
 {
diff --git a/tests/tconst_log2.c b/tests/tconst_log2.c
index 9096cb524..bdfb2db51 100644
--- a/tests/tconst_log2.c
+++ b/tests/tconst_log2.c
@@ -160,7 +160,7 @@ main (int argc, char *argv[])
   tests_start_mpfr ();
 
   p = (argc>1) ? atoi(argv[1]) : 53;
-  rnd = (argc>2) ? (mpfr_rnd_t) atoi(argv[2]) : MPFR_RNDZ;
+  rnd = (argc>2) ? (mpfr_rnd_t) atoi(argv[2]) : MPFR_RNDN;
 
   mpfr_init (x);
 
@@ -180,11 +180,10 @@ main (int argc, char *argv[])
       exit (1);
     }
 
-  if (argc>=2)
+  if (argc >= 2)
     {
       mpfr_set_prec (x, p);
       mpfr_const_log2 (x, rnd);
-      printf ("log(2)=");
       mpfr_out_str (stdout, 10, 0, x, rnd);
       puts ("");
     }
diff --git a/tests/tgeneric.c b/tests/tgeneric.c
index a4189826d..78a852115 100644
--- a/tests/tgeneric.c
+++ b/tests/tgeneric.c
@@ -30,10 +30,16 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
 /* TODO: Add support for type long and extreme integer values, as done
    in tgeneric_ui.c; then tgeneric_ui.c could probably disappear. */
 
+#ifndef ONE_ARG
 #if defined(TWO_ARGS) || defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2) || \
   defined(ULONG_ARG1) || defined(ULONG_ARG2)
 #define TWO_ARGS_ALL
 #endif
+#endif
+
+#if defined(TWO_ARGS_ALL) || defined(ULONG_ARG1) || defined(ULONG_ARG2)
+#define NEED_U
+#endif
 
 #ifndef TEST_RANDOM_POS
 /* For the random function: one number on two is negative. */
@@ -141,7 +147,7 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
 {
   mpfr_prec_t prec, xprec, yprec;
   mpfr_t x, y, z, t, w;
-#if defined(TWO_ARGS_ALL)
+#ifdef NEED_U
   mpfr_t u;
 #endif
 #if defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
@@ -161,7 +167,7 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
   old_emax = mpfr_get_emax ();
 
   mpfr_inits2 (MPFR_PREC_MIN, x, y, z, t, w, (mpfr_ptr) 0);
-#if defined(TWO_ARGS_ALL)
+#ifdef NEED_U
   mpfr_init2 (u, MPFR_PREC_MIN);
 #endif
 
@@ -319,6 +325,8 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
           compare = TEST_FUNCTION (y, x, d, rnd);
           /* d can be infinite due to overflow in mpfr_get_d */
           infinite_input |= DOUBLE_ISINF (d);
+#elif defined(ULONG_ARG1) && defined(ONE_ARG)
+          compare = TEST_FUNCTION (y, i, rnd);
 #elif defined(ULONG_ARG1)
           compare = TEST_FUNCTION (y, i, x, rnd);
 #elif defined(ULONG_ARG2)
@@ -381,6 +389,8 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
                     inexact = TEST_FUNCTION (w, d, x, rnd);
 #elif defined(DOUBLE_ARG2)
                     inexact = TEST_FUNCTION (w, x, d, rnd);
+#elif defined(ULONG_ARG1) && defined(ONE_ARG)
+                    inexact = TEST_FUNCTION (w, i, rnd);
 #elif defined(ULONG_ARG1)
                     inexact = TEST_FUNCTION (w, i, x, rnd);
 #elif defined(ULONG_ARG2)
@@ -400,7 +410,7 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
                                 (mpfr_eexp_t) oemin, (mpfr_eexp_t) e - 1);
                         printf ("x = ");
                         mpfr_dump (x);
-#if defined(TWO_ARGS_ALL)
+#ifdef NEED_U
                         printf ("u = ");
                         mpfr_dump (u);
 #endif
@@ -429,6 +439,8 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
                     inexact = TEST_FUNCTION (w, d, x, rnd);
 #elif defined(DOUBLE_ARG2)
                     inexact = TEST_FUNCTION (w, x, d, rnd);
+#elif defined(ULONG_ARG1) && defined(ONE_ARG)
+                    inexact = TEST_FUNCTION (w, i, rnd);
 #elif defined(ULONG_ARG1)
                     inexact = TEST_FUNCTION (w, i, x, rnd);
 #elif defined(ULONG_ARG2)
@@ -448,7 +460,7 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
                                 (mpfr_eexp_t) e + 1, (mpfr_eexp_t) oemax);
                         printf ("x = ");
                         mpfr_dump (x);
-#if defined(TWO_ARGS_ALL)
+#ifdef NEED_U
                         printf ("u = ");
                         mpfr_dump (u);
 #endif
@@ -492,6 +504,8 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
                 inexact = TEST_FUNCTION (w, d, x, rnd);
 #elif defined(DOUBLE_ARG2)
                 inexact = TEST_FUNCTION (w, x, d, rnd);
+#elif defined(ULONG_ARG1) && defined(ONE_ARG)
+                inexact = TEST_FUNCTION (w, i, rnd);
 #elif defined(ULONG_ARG1)
                 inexact = TEST_FUNCTION (w, i, x, rnd);
 #elif defined(ULONG_ARG2)
@@ -512,7 +526,7 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
                             (mpfr_eexp_t) emin, (mpfr_eexp_t) emax);
                     printf ("x = ");
                     mpfr_dump (x);
-#if defined(TWO_ARGS_ALL)
+#ifdef NEED_U
                     printf ("u = ");
                     mpfr_dump (u);
 #endif
@@ -601,6 +615,8 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
               inexact = TEST_FUNCTION (z, d, x, rnd);
 #elif defined(DOUBLE_ARG2)
               inexact = TEST_FUNCTION (z, x, d, rnd);
+#elif defined(ULONG_ARG1) && defined(ONE_ARG)
+              inexact = TEST_FUNCTION (z, i, rnd);
 #elif defined(ULONG_ARG1)
               inexact = TEST_FUNCTION (z, i, x, rnd);
 #elif defined(ULONG_ARG2)
@@ -615,7 +631,7 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
                   printf ("tgeneric: results differ for "
                           MAKE_STR(TEST_FUNCTION) " on\n  x = ");
                   mpfr_dump (x);
-#if defined(TWO_ARGS_ALL)
+#ifdef NEED_U
                   printf ("  u = ");
                   mpfr_dump (u);
 #endif
@@ -642,7 +658,7 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
                           "\n", mpfr_print_rnd_mode (rnd), compare, inexact);
                   printf ("x = ");
                   mpfr_dump (x);
-#if defined(TWO_ARGS_ALL)
+#ifdef NEED_U
                   printf ("u = ");
                   mpfr_dump (u);
 #endif
@@ -668,7 +684,7 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
                           mpfr_print_rnd_mode (rnd));
                   printf ("x = ");
                   mpfr_dump (x);
-#if defined(TWO_ARGS_ALL)
+#ifdef NEED_U
                   printf ("u = ");
                   mpfr_dump (u);
 #endif
@@ -691,7 +707,7 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
 #endif
 
   mpfr_clears (x, y, z, t, w, (mpfr_ptr) 0);
-#if defined(TWO_ARGS_ALL)
+#ifdef NEED_U
   mpfr_clear (u);
 #endif
 }
@@ -704,6 +720,7 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
 #undef RAND_FUNCTION
 #undef TWO_ARGS
 #undef TWO_ARGS_ALL
+#undef NEED_U
 #undef DOUBLE_ARG1
 #undef DOUBLE_ARG2
 #undef ULONG_ARG1
diff --git a/tests/tlog.c b/tests/tlog.c
index fc1e3fba9..21d2933c6 100644
--- a/tests/tlog.c
+++ b/tests/tlog.c
@@ -67,18 +67,16 @@ check2 (const char *as, mpfr_rnd_t rnd_mode, const char *res1s)
 }
 
 static void
-check3 (double d, unsigned long prec, mpfr_rnd_t rnd)
+check3 (char *s, unsigned long prec, mpfr_rnd_t rnd)
 {
   mpfr_t x, y;
 
   mpfr_init2 (x, prec);
   mpfr_init2 (y, prec);
-  mpfr_set_d (x, d, rnd);
+  mpfr_set_str (x, s, 10, rnd);
   test_log (y, x, rnd);
   mpfr_out_str (stdout, 10, 0, y, rnd);
   puts ("");
-  mpfr_print_binary (y);
-  puts ("");
   mpfr_clear (x);
   mpfr_clear (y);
 }
@@ -274,9 +272,10 @@ main (int argc, char *argv[])
 {
   tests_start_mpfr ();
 
-  if (argc==4)
+  if (argc == 3 || argc == 4)
     {   /* tlog x prec rnd */
-      check3 (atof(argv[1]), atoi(argv[2]), (mpfr_rnd_t) atoi(argv[3]));
+      check3 (argv[1], strtoul (argv[2], NULL, 10),
+              (argc == 4) ? atoi (argv[3]) : MPFR_RNDN);
       goto done;
     }
 
diff --git a/tests/tlog_ui.c b/tests/tlog_ui.c
index 50d4f3658..ca145c27d 100644
--- a/tests/tlog_ui.c
+++ b/tests/tlog_ui.c
@@ -22,6 +22,33 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
 
 #include "mpfr-test.h"
 
+static void
+compare_with_log (unsigned long n, mpfr_prec_t p)
+{
+  mpfr_t x, y;
+  int inex1, inex2;
+  mpfr_flags_t flags1;
+
+  mpfr_init2 (x, sizeof (unsigned long) * CHAR_BIT);
+  mpfr_init2 (y, p);
+  inex1 = mpfr_set_ui (x, n, MPFR_RNDN);
+  MPFR_ASSERTN(inex1 == 0);
+  inex1 = mpfr_log (y, x, MPFR_RNDN);
+  flags1 = __gmpfr_flags;
+  mpfr_set_prec (x, p);
+  inex2 = mpfr_log_ui (x, n, MPFR_RNDN);
+  MPFR_ASSERTN(inex1 == inex2);
+  MPFR_ASSERTN(flags1 == __gmpfr_flags);
+  MPFR_ASSERTN(mpfr_equal_p (x, y));
+  mpfr_clears (x, y, (mpfr_ptr) 0);
+}
+
+#define TEST_FUNCTION mpfr_log_ui
+#define ONE_ARG
+#define ULONG_ARG1
+#define RAND_FUNCTION(x) mpfr_set_ui (x, randlimb (), MPFR_RNDN)
+#include "tgeneric.c"
+
 #define TEST_FUNCTION mpfr_log_ui
 
 int
@@ -45,8 +72,8 @@ main (int argc, char *argv[])
 
   if (argc >= 3) /* tlog_ui n prec [rnd] */
     {
-      mpfr_set_prec (x, atoi (argv[2]));
-      TEST_FUNCTION (x, atoi (argv[1]),
+      mpfr_set_prec (x, strtoul (argv[2], NULL, 10));
+      TEST_FUNCTION (x, strtoul (argv[1], NULL, 10),
                     argc > 3 ? (mpfr_rnd_t) atoi (argv[3]) : MPFR_RNDN);
       mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
       printf ("\n");
@@ -161,6 +188,12 @@ main (int argc, char *argv[])
           }
     }
 
+
+  test_generic (MPFR_PREC_MIN, 1000, 1);
+
+  for (n = 1; n < ULONG_MAX / 3; n *= 3)
+    compare_with_log (n, 10);
+
  clear_and_exit:
   mpfr_clears (x, y, z, t, v, (mpfr_ptr) 0);
 
diff --git a/tests/tsub.c b/tests/tsub.c
index a37399d81..7ad321e5a 100644
--- a/tests/tsub.c
+++ b/tests/tsub.c
@@ -629,6 +629,158 @@ check_rounding (void)
     }
 }
 
+/* Check a = b - c, where the significand of b has all 1's, c is small
+   compared to b, and PREC(a) = PREC(b) - 1. Thus b is a midpoint for
+   the precision of the result a. The test is done with the extended
+   exponent range and with some reduced exponent range. Two choices
+   are made for the exponent of b: the maximum exponent - 1 (similar
+   to some normal case) and the maximum exponent (overflow case or
+   near overflow case, depending on the rounding mode).
+   This test is useful to trigger a bug in r10382: Since c is small,
+   the computation in sub1.c was done by first rounding b in the
+   precision of a, then correcting the result if b was a breakpoint
+   for this precision (exactly representable number for the directed
+   rounding modes, or midpoint for the round-to-nearest mode). The
+   problem was that for a midpoint in the round-to-nearest mode, the
+   rounding of b gave a spurious overflow; not only the overflow flag
+   was incorrect, but the result could not be corrected, since due to
+   this overflow, the "even rounding" information was lost.
+   In the case of reduced exponent range, an additional test is done
+   for consistency checks: the subtraction is done in the extended
+   exponent range (no overflow), then the result is converted to the
+   initial exponent range with mpfr_check_range. */
+static void
+check_max_almosteven (void)
+{
+  mpfr_exp_t old_emin, old_emax;
+  mpfr_exp_t emin[2] = { MPFR_EMIN_MIN, -1000 };
+  mpfr_exp_t emax[2] = { MPFR_EMAX_MAX, 1000 };
+  int i;
+
+  old_emin = mpfr_get_emin ();
+  old_emax = mpfr_get_emax ();
+
+  for (i = 0; i < 2; i++)
+    {
+      mpfr_t a1, a2, b, c;
+      mpfr_prec_t p;
+      int neg, j, rnd;
+
+      set_emin (emin[i]);
+      set_emax (emax[i]);
+
+      p = MPFR_PREC_MIN + randlimb () % 70;
+      mpfr_init2 (a1, p);
+      mpfr_init2 (a2, p);
+      mpfr_init2 (b, p+1);
+      mpfr_init2 (c, MPFR_PREC_MIN);
+
+      mpfr_setmax (b, 0);
+      mpfr_set_ui (c, 1, MPFR_RNDN);
+
+      for (neg = 0; neg < 2; neg++)
+        {
+          for (j = 1; j >= 0; j--)
+            {
+              mpfr_set_exp (b, __gmpfr_emax - j);
+              RND_LOOP (rnd)
+                {
+                  mpfr_flags_t flags1, flags2;
+                  int inex1, inex2;
+
+                  /* Expected result. */
+                  flags1 = MPFR_FLAGS_INEXACT;
+                  if (rnd == MPFR_RNDN || MPFR_IS_LIKE_RNDZ (rnd, neg))
+                    {
+                      inex1 = neg ? 1 : -1;
+                      mpfr_setmax (a1, __gmpfr_emax - j);
+                    }
+                  else
+                    {
+                      inex1 = neg ? -1 : 1;
+                      if (j == 0)
+                        {
+                          flags1 |= MPFR_FLAGS_OVERFLOW;
+                          mpfr_set_inf (a1, 1);
+                        }
+                      else
+                        {
+                          mpfr_setmin (a1, __gmpfr_emax);
+                        }
+                    }
+                  MPFR_SET_SIGN (a1, neg ? -1 : 1);
+
+                  /* Computed result. */
+                  mpfr_clear_flags ();
+                  inex2 = mpfr_sub (a2, b, c, (mpfr_rnd_t) rnd);
+                  flags2 = __gmpfr_flags;
+
+                  if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
+                         mpfr_equal_p (a1, a2)))
+                    {
+                      printf ("Error 1 in check_max_almosteven for %s,"
+                              " i = %d, j = %d, neg = %d\n",
+                              mpfr_print_rnd_mode ((mpfr_rnd_t) rnd),
+                              i, j, neg);
+                      printf ("     b = ");
+                      mpfr_dump (b);
+                      printf ("Expected ");
+                      mpfr_dump (a1);
+                      printf ("  with inex = %d, flags =", inex1);
+                      flags_out (flags1);
+                      printf ("Got      ");
+                      mpfr_dump (a2);
+                      printf ("  with inex = %d, flags =", inex2);
+                      flags_out (flags2);
+                      exit (1);
+                    }
+
+                  if (i == 0)
+                    break;
+
+                  /* Additional test for the reduced exponent range. */
+                  mpfr_clear_flags ();
+                  set_emin (MPFR_EMIN_MIN);
+                  set_emax (MPFR_EMAX_MAX);
+                  inex2 = mpfr_sub (a2, b, c, (mpfr_rnd_t) rnd);
+                  set_emin (emin[i]);
+                  set_emax (emax[i]);
+                  inex2 = mpfr_check_range (a2, inex2, (mpfr_rnd_t) rnd);
+                  flags2 = __gmpfr_flags;
+
+                  if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
+                         mpfr_equal_p (a1, a2)))
+                    {
+                      printf ("Error 2 in check_max_almosteven for %s,"
+                              " i = %d, j = %d, neg = %d\n",
+                              mpfr_print_rnd_mode ((mpfr_rnd_t) rnd),
+                              i, j, neg);
+                      printf ("     b = ");
+                      mpfr_dump (b);
+                      printf ("Expected ");
+                      mpfr_dump (a1);
+                      printf ("  with inex = %d, flags =", inex1);
+                      flags_out (flags1);
+                      printf ("Got      ");
+                      mpfr_dump (a2);
+                      printf ("  with inex = %d, flags =", inex2);
+                      flags_out (flags2);
+                      exit (1);
+                    }
+                }
+            }  /* j */
+
+          mpfr_neg (b, b, MPFR_RNDN);
+          mpfr_neg (c, c, MPFR_RNDN);
+        }  /* neg */
+
+      mpfr_clears (a1, a2, b, c, (mpfr_ptr) 0);
+    }  /* i */
+
+  set_emin (old_emin);
+  set_emax (old_emax);
+}
+
 #define TEST_FUNCTION test_sub
 #define TWO_ARGS
 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS)
@@ -646,6 +798,7 @@ main (void)
   check_rounding ();
   check_diverse ();
   check_inexact ();
+  check_max_almosteven ();
   bug_ddefour ();
   for (p=2; p<200; p++)
     for (i=0; i<50; i++)
diff --git a/tools/nightly-test b/tools/nightly-test
index 8134f9cd0..45882535d 100755
--- a/tools/nightly-test
+++ b/tools/nightly-test
@@ -15,7 +15,7 @@ mkdir "$DIR"
 test ! -h "$DIR"
 svn checkout "svn://scm.gforge.inria.fr/svn/mpfr/${BRANCH:-trunk}" "$DIR"
 cd "$DIR"
-autoreconf -i
+./autogen.sh
 # -Wmissing-prototypes is not valid for C++
 # add -Wno-sign-conversion because of a bug in g++
 # -pedantic-errors (instead of just -pedantic) below allows us to