Replaced the 4.0 branch by a copy of the trunk (in order to keep the

history of the latest changes of the trunk); this replacement is due
to a limitation of Subversion where a merge is regarded as a single
commit, and is OK here since the 4.0 branch is new.

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

16 files changed, 391 insertions, 127 deletions

diff --git a/NEWS b/NEWS
index 21cff1f39..76a03b3d0 100644
--- a/NEWS
+++ b/NEWS
@@ -110,7 +110,6 @@ Changes from versions 3.1.* to version 4.0.0:
   in the computation of Bernoulli numbers (used in mpfr_gamma, mpfr_li2,
   mpfr_digamma, mpfr_lngamma and mpfr_lgamma), in mpfr_div, in mpfr_fma
   and mpfr_fms.
-- Test coverage: 96.3% lines of code.
 - Bug fixes. In particular: a speed improvement when the --enable-assert
   or --enable-assert=full configure option is used with GCC; mpfr_get_str
   now sets the NaN flag on NaN input and the inexact flag when the conversion
diff --git a/doc/FAQ.html b/doc/FAQ.html
index ae9f29fae..92530c190 100644
--- a/doc/FAQ.html
+++ b/doc/FAQ.html
@@ -338,7 +338,7 @@ installation using <cite>autoconf</cite> or <cite>pkg-config</cite>?</dt>
 <dd><p>The <cite><acronym>MPFR</acronym></cite> team does not currently
 recommend any <cite>autoconf</cite> code, but a section will later
 be added to the <cite><acronym>MPFR</acronym></cite> manual.
-Limited <cite>pkg-config</cite> support has been added for
+Limited <cite>pkg-config</cite> support will be added for
 <cite><acronym>MPFR</acronym></cite> 4.0.0; example:</p>
 <pre style="margin-left: 2em">cc myprogram.c $(pkg-config --cflags --libs mpfr)</pre></dd>
 
diff --git a/doc/algorithms.tex b/doc/algorithms.tex
index d5311d1c3..3010a2738 100644
--- a/doc/algorithms.tex
+++ b/doc/algorithms.tex
@@ -1854,37 +1854,47 @@ that $x \geq 0$.
 
 We use Higham's notation, with $\theta_i$ denoting variables such that
 $|\theta_i| \leq 2^{-p}$.
-Firstly, $u$ is exact. Then $v = e^{2x} (1+\theta_1)$ and
+Firstly, $u$ is exact, assuming $x$ is exact with precision $p$.
+Then $v = e^{2x} (1+\theta_1)$ and
 $w = (e^{2x}+1) (1+\theta_2)^2$.
 The error on $r$ is bounded by $\frac{1}{2} \ulp(v) + \frac{1}{2} \ulp(r)$.
-Assume $\ulp(v) = 2^e \ulp(r)$, with $e \geq 0$;
-then the error on $r$ is bounded by $\frac{1}{2} (2^e+1) \ulp(r)$.
-We can thus write $r = (e^{2x}-1) (1+\theta_3)^{2^e+1}$,
-and then $s = \tanh(x) (1+\theta_4)^{2^e+4}$.
+Assume $\ulp(v) = 2^k \ulp(r)$, with $k \geq 0$;
+then the error on $r$ is bounded by $\frac{1}{2} (2^k+1) \ulp(r)$.
+We can thus write $r = (e^{2x}-1) (1+\theta_3)^{2^k+1}$,
+and then $s = \tanh(x) \cdot (1+\theta_4)^{2^k+4}$.
 
 \begin{lemma}
-For $|x| \leq 1/2$, and $|y| \leq |x|^{-1/2}$, we have:
-\[ |(1+x)^y-1| \leq 2 |y| x. \]
+For $0 < x \leq 1/2$ and $0 < y \leq x^{-1/2}$, we have:
+\[ 0 < (1+x)^y - 1 \leq 1.4 \cdot y \cdot x. \]
 \end{lemma}
 \begin{proof}
-We have $(1+x)^y = e^{y \log (1+x)}$,
-with $|y \log (1+x)| \leq |x|^{-1/2} |\log (1+x)|$.
-The function $|x|^{-1/2} \log (1+x)$ is increasing on $[-1/2,1/2]$, and
-takes as values $\approx -0.490$ in $x=-1/2$ and $\approx 0.286$ in $x=1/2$,
-thus is bounded in absolute value by $1/2$.
-This yields $|y \log (1+x)| \leq 1/2$.
-Now it is easy to see that for $|t| \leq 1/2$, we have
-$|e^t-1| \leq 1.3 |t|$.
-Thus $|(1+x)^y-1| \leq 1.3 |y| |\log (1+x)|$.
-The result follows from $|\log (1+x)| \leq 1.4 |x|$ for $|x| \leq 1/2$,
-and $1.3 \times 1.4 \leq 2$.
+We have $(1+x)^y = e^{y \cdot \log(1+x)}$, with
+$0 < y \cdot \log(1+x) \leq x^{-1/2} \cdot \log(1+x)$.
+The function $x^{-1/2} \cdot \log(1+x)$ is increasing on $]0,1/2]$,
+and reaches $\approx 0.573$ for $x = 1/2$.
+Thus $0 < y \cdot \log(1+x) < 0.574$.
+Now it is easy to see that for $0 < t < 0.574$, we have
+$|e^t-1| \leq 1.4\,t$.
+Thus $0 < (1+x)^y - 1 \leq 1.4 \cdot y \cdot \log(1+x)$.
+The result follows from $\log(1+x) \leq x$ for
+$0 < x \leq 1/2$.
 \end{proof}
 
-Applying the above lemma for $x=\theta_4$ and $y=2^e+4$,
-assuming $2^e+4 \leq 2^{p/2}$,
-we get $s = \tanh(x) [1 + 2(2^e+4)\theta_5]$.
-Since $2^e+4 \leq 2^{{\rm max}(3,e+1)}$,
-the relative error on $s$ is thus bounded by $2^{{\rm max}(4,e+2)-p}$.
+First, note that for $t > 0$ and $q \geq 1$, one has
+$|(1-t)^q - 1| \leq |(1+t)^q - 1|$ due to the triangle inequality on the
+development. Thus $|(1+\theta_4)^{2^k+4} - 1| \leq |(1+2^{-p})^{2^k+4} - 1|$.
+Then one can apply the above lemma for $x=2^{-p}$ and $y=2^k+4$,
+assuming $2^k+4 \leq 2^{p/2}$.
+We get $|(1+\theta_4)^{2^k+4} - 1| \leq |(1+2^{-p})^{2^k+4} - 1|
+\leq 1.4 (2^k+4) 2^{-p}$, and thus
+we can write $s = \tanh(x) [1 + 1.4 (2^k+4)\theta_5]$ with
+$|\theta_5| \leq 2^{-p}$.
+Since $2^k+4 \leq 2^{{\rm max}(3,k+1)}$,
+the relative error on $s$ is thus bounded by $2^{{\rm max}(4,k+2)-p}$.
+
+The condition $2^k+4 \leq 2^{p/2}$ is checked in the code with
+${\rm max}(3,k+1) \leq \lfloor p/2 \rfloor$, so that:
+\[2^k+4 \leq 2^{{\rm max}(3,k+1)} \leq 2^{\lfloor p/2 \rfloor} \leq 2^{p/2}.\]
 
 \subsection{The inverse hyperbolic tangent function}
 
diff --git a/doc/mpfr.texi b/doc/mpfr.texi
index 3dccf760b..01da35cec 100644
--- a/doc/mpfr.texi
+++ b/doc/mpfr.texi
@@ -15,9 +15,9 @@
 @c as the "info" utility cannot currently handle them.
 @c   https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=212549
 
-@c texinfo.tex 2017-12-01.13 or later is needed to support @var{} in exponent.
-@c See thread at
-@c   http://lists.gnu.org/archive/html/bug-texinfo/2017-11/msg00019.html
+@c Warning! If a macro is used with @iftex, it must also be defined for
+@c info at least (e.g. with @ifnottex) in order to avoid a failure with
+@c makeinfo 5.2 (there is no such problem with makeinfo 6.5).
 
 @copying
 This manual documents how to install and use the Multiple Precision
@@ -107,6 +107,30 @@ floating-point arithmetic, version @value{VERSION}.
 @end menu
 
 
+@c Texinfo currently doesn't support @var in script font size:
+@c error "@scriptfont 5 is undefined". See thread at
+@c   http://lists.gnu.org/archive/html/bug-texinfo/2017-11/msg00019.html
+@c   http://lists.gnu.org/archive/html/bug-texinfo/2017-12/msg00001.html
+@c (unfortunately the antispam there garbles important Texinfo code).
+@c
+@c Use the @svar macro, defined below, instead. It currently disables
+@c the use of @var, with the consequence that this appears in italic
+@c instead of slanted. Alternatively, the following could be used for
+@c the macro definition: @hbox{@switchtolllsize@var{\V\}}
+@c But with this workaround, the font is larger than it should be, making
+@c some cases ugly, such as in: {\log 2 \over \log @svar{b}}
+@iftex
+@macro svar {V}
+\V\
+@end macro
+@end iftex
+@ifnottex
+@macro svar {V}
+@var{\V\}
+@end macro
+@end ifnottex
+
+
 @c  @m{T,N} is $T$ in tex or @math{N} otherwise.  This is an easy way to give
 @c  different forms for math in tex and info.  Commas in N or T don't work,
 @c  but @C{} can be used instead.  \, works in info but not in tex.
@@ -1690,7 +1714,7 @@ minimal precision @math{m} depending only on @var{p} = PREC(@var{op}) and
 m = 1 + ceil(@var{p}*log(2)/log(@var{b})),
 @end ifnottex
 @tex
-$m = 1 + \left\lceil @var{p} {\log 2 \over \log @var{b}} \right\rceil$,
+$m = 1 + \left\lceil @var{p} {\log 2 \over \log @svar{b}} \right\rceil$,
 @end tex
 with @var{p} replaced by @var{p}@minus{}1 if @var{b} is a power of 2,
 but in some very rare cases, it might be @math{m+1}
@@ -1965,7 +1989,7 @@ or both arguments are NaN@. But only floating-point numbers can be compared
 
 @deftypefun int mpfr_cmp_ui_2exp (mpfr_t @var{op1}, unsigned long int @var{op2}, mpfr_exp_t @var{e})
 @deftypefunx int mpfr_cmp_si_2exp (mpfr_t @var{op1}, long int @var{op2}, mpfr_exp_t @var{e})
-Compare @var{op1} and @m{@var{op2} \times 2^{@var{e}}, @var{op2} multiplied by two to
+Compare @var{op1} and @m{@var{op2} \times 2^{@svar{e}}, @var{op2} multiplied by two to
 the power @var{e}}. Similar as above.
 @end deftypefun
 
@@ -2205,7 +2229,7 @@ at @minus{}@var{op} (formula 5.1.1 from the same reference).
 @deftypefun int mpfr_li2 (mpfr_t @var{rop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
 Set @var{rop} to real part of the dilogarithm of @var{op}, rounded in the
 direction @var{rnd}. MPFR defines the dilogarithm function as
-@m{-\int_{t=0}^{@var{op}} \log(1-t)/t\ dt,the integral of -log(1-t)/t from 0
+@m{-\int_{t=0}^{@svar{op}} \log(1-t)/t\ dt,the integral of -log(1-t)/t from 0
 to @var{op}}.
 @end deftypefun
 
diff --git a/doc/texinfo.tex b/doc/texinfo.tex
index 5d05ba5ca..42bb98335 100644
--- a/doc/texinfo.tex
+++ b/doc/texinfo.tex
@@ -3,7 +3,7 @@
 % Load plain if necessary, i.e., if running under initex.
 \expandafter\ifx\csname fmtname\endcsname\relax\input plain\fi
 %
-\def\texinfoversion{2017-12-01.13}
+\def\texinfoversion{2017-12-18.20}
 %
 % Copyright 1985, 1986, 1988, 1990, 1991, 1992, 1993, 1994, 1995,
 % 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006,
@@ -2503,25 +2503,13 @@ end
 
 
 % In order for the font changes to affect most math symbols and letters,
-% we have to define the \textfont of the standard families.
+% we have to define the \textfont of the standard families.  We don't
+% bother to reset \scriptfont and \scriptscriptfont; awaiting user need.
 %
 \def\resetmathfonts{%
   \textfont0=\rmfont \textfont1=\ifont \textfont2=\syfont
   \textfont\itfam=\itfont \textfont\slfam=\slfont \textfont\bffam=\bffont
   \textfont\ttfam=\ttfont \textfont\sffam=\sffont
-  %
-  \scriptfont0=\rmfont \scriptfont1=\ifont \scriptfont2=\syfont
-  \scriptfont\itfam=\itfont \scriptfont\slfam=\slfont \scriptfont\bffam=\bffont
-  \scriptfont\ttfam=\ttfont \scriptfont\sffam=\sffont
-  %
-  \scriptscriptfont0=\rmfont
-  \scriptscriptfont1=\ifont
-  \scriptscriptfont2=\syfont
-  \scriptscriptfont\itfam=\itfont
-  \scriptscriptfont\slfam=\slfont
-  \scriptscriptfont\bffam=\bffont
-  \scriptscriptfont\ttfam=\ttfont
-  \scriptscriptfont\sffam=\sffont
 }
 
 %
diff --git a/src/div.c b/src/div.c
index 3ef62b6d5..7d208621f 100644
--- a/src/div.c
+++ b/src/div.c
@@ -71,7 +71,8 @@ mpfr_div2_approx (mpfr_limb_ptr Q1, mpfr_limb_ptr Q0,
 
   /* we ignore yy below, but first increment r0, to ensure we get a lower
      approximation of the remainder */
-  r0 += (yy != 0);
+  r0 += yy != 0;
+  r1 += r0 == 0 && yy != 0;
   r0 = u0 - r0;
   r1 = u1 - r1 - (r0 > u0);
 
diff --git a/src/exp2.c b/src/exp2.c
index 0badfaf3b..8aa133224 100644
--- a/src/exp2.c
+++ b/src/exp2.c
@@ -29,7 +29,7 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
 int
 mpfr_exp2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
 {
-  int inexact;
+  int inexact, inex2;
   long xint;
   mpfr_t xfrac;
   MPFR_SAVE_EXPO_DECL (expo);
@@ -62,24 +62,22 @@ mpfr_exp2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
         }
     }
 
-  /* since the smallest representable non-zero float is 1/2*2^__gmpfr_emin,
-     if x < __gmpfr_emin - 1, the result is either 1/2*2^__gmpfr_emin or 0 */
-  MPFR_ASSERTN (MPFR_EMIN_MIN >= LONG_MIN + 2);
-  if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emin - 1) < 0))
-    {
-      mpfr_rnd_t rnd2 = rnd_mode;
-      /* in round to nearest mode, round to zero when x <= __gmpfr_emin-2 */
-      if (rnd_mode == MPFR_RNDN &&
-          mpfr_cmp_si_2exp (x, __gmpfr_emin - 2, 0) <= 0)
-        rnd2 = MPFR_RNDZ;
-      return mpfr_underflow (y, rnd2, 1);
-    }
+  /* Since the smallest representable non-zero float is 1/2 * 2^emin,
+     if x <= emin - 2, the result is either 1/2 * 2^emin or 0.
+     Warning, for emin - 2 < x < emin - 1, we cannot conclude, since 2^x
+     might round to 2^(emin - 1) for rounding away or to nearest, and there
+     might be no underflow, since we consider underflow "after rounding". */
+
+  MPFR_STAT_STATIC_ASSERT (MPFR_EMIN_MIN >= LONG_MIN + 2);
+  if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emin - 2) <= 0))
+    return mpfr_underflow (y, rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode, 1);
 
-  MPFR_ASSERTN (MPFR_EMAX_MAX <= LONG_MAX);
+  MPFR_STAT_STATIC_ASSERT (MPFR_EMAX_MAX <= LONG_MAX);
   if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emax) >= 0))
     return mpfr_overflow (y, rnd_mode, 1);
 
-  /* We now know that emin - 1 <= x < emax. */
+  /* We now know that emin - 2 < x < emax. Note that an underflow or
+     overflow is still possible (we have eliminated only easy cases). */
 
   MPFR_SAVE_EXPO_MARK (expo);
 
@@ -91,10 +89,13 @@ mpfr_exp2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
 
   xint = mpfr_get_si (x, MPFR_RNDZ);
   mpfr_init2 (xfrac, MPFR_PREC (x));
-  mpfr_sub_si (xfrac, x, xint, MPFR_RNDN); /* exact */
+  MPFR_DBGRES (inexact = mpfr_sub_si (xfrac, x, xint, MPFR_RNDN));
+  MPFR_ASSERTD (inexact == 0);
 
   if (MPFR_IS_ZERO (xfrac))
     {
+      /* Here, emin - 1 <= x <= emax - 1, so that an underflow or overflow
+         will not be possible. */
       mpfr_set_ui (y, 1, MPFR_RNDN);
       inexact = 0;
     }
@@ -141,11 +142,27 @@ mpfr_exp2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
     }
 
   mpfr_clear (xfrac);
-  MPFR_CLEAR_FLAGS ();
-  mpfr_mul_2si (y, y, xint, MPFR_RNDN); /* exact or overflow */
-  /* Note: We can have an overflow only when t was rounded up to 2. */
-  MPFR_ASSERTD (MPFR_IS_PURE_FP (y) || inexact > 0);
-  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
+
+  if (MPFR_UNLIKELY (rnd_mode == MPFR_RNDN && xint == __gmpfr_emin - 1 &&
+                     MPFR_GET_EXP (y) == 0 && mpfr_powerof2_raw (y)))
+    {
+      /* y was rounded down to 1/2 and the rounded value with an unbounded
+         exponent range would be 2^(emin-2), i.e. the midpoint between 0
+         and the smallest positive FP number. This is a double rounding
+         problem: we should not round to 0, but to (1/2) * 2^emin. */
+      MPFR_SET_EXP (y, __gmpfr_emin);
+      inexact = 1;
+      MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_UNDERFLOW);
+    }
+  else
+    {
+      MPFR_CLEAR_FLAGS ();
+      inex2 = mpfr_mul_2si (y, y, xint, rnd_mode);
+      if (inex2 != 0)  /* underflow or overflow */
+        inexact = inex2;
+      MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
+    }
+
   MPFR_SAVE_EXPO_FREE (expo);
   return mpfr_check_range (y, inexact, rnd_mode);
 }
diff --git a/src/fma.c b/src/fma.c
index 958176035..5acbec69f 100644
--- a/src/fma.c
+++ b/src/fma.c
@@ -103,6 +103,7 @@ mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z,
   mpfr_t u;
   mp_size_t n;
   mpfr_exp_t e;
+  mpfr_prec_t precx, precy;
   MPFR_SAVE_EXPO_DECL (expo);
   MPFR_GROUP_DECL(group);
 
@@ -121,22 +122,27 @@ mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z,
 
   e = MPFR_GET_EXP (x) + MPFR_GET_EXP (y);
 
+  precx = MPFR_PREC (x);
+  precy = MPFR_PREC (y);
+
   /* First deal with special case where prec(x) = prec(y) and x*y does
      not overflow nor underflow. Do it only for small sizes since for large
      sizes x*y is faster using Mulders' algorithm (as a rule of thumb,
      we assume mpn_mul_n is faster up to 4*MPFR_MUL_THRESHOLD).
      Since |EXP(x)|, |EXP(y)| < 2^(k-2) on a k-bit computer,
      |EXP(x)+EXP(y)| < 2^(k-1), thus cannot overflow nor underflow. */
-  if (MPFR_PREC(x) == MPFR_PREC(y) && e <= __gmpfr_emax && e > __gmpfr_emin)
+  if (precx == precy && e <= __gmpfr_emax && e > __gmpfr_emin)
     {
-      if (MPFR_PREC(x) < GMP_NUMB_BITS && MPFR_PREC(z) == MPFR_PREC(x))
+      if (precx < GMP_NUMB_BITS &&
+          MPFR_PREC(z) == precx &&
+          MPFR_PREC(s) == precx)
         {
           mp_limb_t umant[2], zmant[2];
           mpfr_t zz;
           int inex;
 
           umul_ppmm (umant[1], umant[0], MPFR_MANT(x)[0], MPFR_MANT(y)[0]);
-          MPFR_PREC(u) = MPFR_PREC(zz) = 2 * MPFR_PREC(x);
+          MPFR_PREC(u) = MPFR_PREC(zz) = 2 * precx;
           MPFR_MANT(u) = umant;
           MPFR_MANT(zz) = zmant;
           MPFR_SIGN(u) = MPFR_MULT_SIGN( MPFR_SIGN(x) , MPFR_SIGN(y) );
@@ -162,7 +168,8 @@ mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z,
               zmant[0] = MPFR_LIMB_ZERO;
               if ((umant[1] & MPFR_LIMB_HIGHBIT) == 0)
                 {
-                  umant[1] = (umant[1] << 1) | (umant[0] >> (GMP_NUMB_BITS - 1));
+                  umant[1] = (umant[1] << 1) |
+                    (umant[0] >> (GMP_NUMB_BITS - 1));
                   umant[0] = umant[0] << 1;
                   MPFR_EXP(u) = e - 1;
                 }
@@ -170,6 +177,8 @@ mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z,
                 MPFR_EXP(u) = e;
             }
           inex = mpfr_add (u, u, zz, rnd_mode);
+          /* mpfr_set_1_2 requires PREC(u) = 2*PREC(s),
+             thus we need PREC(s) = PREC(x) = PREC(y) = PREC(z) */
           return mpfr_set_1_2 (s, u, rnd_mode, inex);
         }
       else if ((n = MPFR_LIMB_SIZE(x)) <= 4 * MPFR_MUL_THRESHOLD)
@@ -201,7 +210,8 @@ mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z,
 
   /* If we take prec(u) >= prec(x) + prec(y), the product u <- x*y
      is exact, except in case of overflow or underflow. */
-  MPFR_GROUP_INIT_1 (group, MPFR_PREC(x) + MPFR_PREC(y), u);
+  MPFR_ASSERTN (precx + precy <= MPFR_PREC_MAX);
+  MPFR_GROUP_INIT_1 (group, precx + precy, u);
   MPFR_SAVE_EXPO_MARK (expo);
 
   if (MPFR_UNLIKELY (mpfr_mul (u, x, y, MPFR_RNDN)))
@@ -328,13 +338,13 @@ mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z,
               MPFR_BLOCK (flags,
                           if (MPFR_GET_EXP (x) < MPFR_GET_EXP (y))
                             {
-                              mpfr_init2 (scaled_v, MPFR_PREC (x));
+                              mpfr_init2 (scaled_v, precx);
                               mpfr_mul_2ui (scaled_v, x, scale, MPFR_RNDN);
                               mpfr_mul (u, scaled_v, y, MPFR_RNDN);
                             }
                           else
                             {
-                              mpfr_init2 (scaled_v, MPFR_PREC (y));
+                              mpfr_init2 (scaled_v, precy);
                               mpfr_mul_2ui (scaled_v, y, scale, MPFR_RNDN);
                               mpfr_mul (u, x, scaled_v, MPFR_RNDN);
                             });
diff --git a/src/mul.c b/src/mul.c
index c22a49222..54f8d8170 100644
--- a/src/mul.c
+++ b/src/mul.c
@@ -847,7 +847,7 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
   else
     /* Mulders' mulhigh. This code can also be used via mpfr_sqr,
        hence the tests b != c. */
-    if (MPFR_UNLIKELY (bn > (threshold = b != c ?
+    if (MPFR_UNLIKELY (cn > (threshold = b != c ?
                              MPFR_MUL_THRESHOLD : MPFR_SQR_THRESHOLD)))
       {
         mp_limb_t *bp, *cp;
@@ -981,23 +981,15 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
         /* now the approximation is in tmp[tn-n]...tmp[tn-1] */
         MPFR_ASSERTD (MPFR_LIMB_MSB (tmp[tn-1]) != 0);
 
+        /* for RNDF, we simply use RNDZ, since anyway here we multiply numbers
+           with large precisions, thus the overhead of RNDZ is small */
+        if (rnd_mode == MPFR_RNDF)
+          rnd_mode = MPFR_RNDZ;
+
         /* if the most significant bit b1 is zero, we have only p-1 correct
            bits */
-        if (rnd_mode == MPFR_RNDF && p + b1 - 1 >= aq
-            && ((aq % GMP_NUMB_BITS) != 0))
-          {
-            /* Warning: while we know the error is bounded by
-               ufp(tmp)/2^(p+b1-1), and we want aq correct bits,
-               simply truncating tmp might give a result which is
-               1 ulp too small, in case tmp ends with 111...111.
-               But if we add one to tmp[tn-n], and then truncate, we
-               cannot exceed what would be obtained by RNDU,
-               since tmp has more bits than a */
-            mpn_add_1 (tmp + tn - n, tmp + tn - n, n, MPFR_LIMB_ONE);
-            /* we can round */
-          }
-        else if (MPFR_UNLIKELY (!mpfr_round_p (tmp, tn, p + b1 - 1,
-                                               aq + (rnd_mode == MPFR_RNDN))))
+        if (MPFR_UNLIKELY (!mpfr_round_p (tmp, tn, p + b1 - 1,
+                                          aq + (rnd_mode == MPFR_RNDN))))
           {
             tmp -= k - tn; /* tmp may have changed, FIX IT!!!!! */
             goto full_multiply;
diff --git a/src/tanh.c b/src/tanh.c
index 06adf7874..5459ecfa4 100644
--- a/src/tanh.c
+++ b/src/tanh.c
@@ -98,50 +98,62 @@ mpfr_tanh (mpfr_ptr y, mpfr_srcptr xt , mpfr_rnd_t rnd_mode)
     if (MPFR_GET_EXP (x) < 0)
       Nt += -MPFR_GET_EXP (x);
 
+    /* The error analysis in algorithms.tex assumes that x is exactly
+       representable with working precision Nt.
+       FIXME: adapt the error analysis for the case Nt < PREC(x). */
+    if (Nt < MPFR_PREC(x))
+      Nt = MPFR_PREC(x);
+
     /* initialize of intermediary variable */
     MPFR_GROUP_INIT_2 (group, Nt, t, te);
 
     MPFR_ZIV_INIT (loop, Nt);
-    for (;;) {
-      /* tanh = (exp(2x)-1)/(exp(2x)+1) */
-      mpfr_mul_2ui (te, x, 1, MPFR_RNDN);  /* 2x */
-      /* since x > 0, we can only have an overflow */
-      mpfr_exp (te, te, MPFR_RNDN);        /* exp(2x) */
-      if (MPFR_UNLIKELY (MPFR_IS_INF (te))) {
-      set_one:
-        inexact = MPFR_FROM_SIGN_TO_INT (sign);
-        mpfr_set4 (y, __gmpfr_one, MPFR_RNDN, sign);
-        if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG_SIGN (sign)))
+    for (;;)
+      {
+        /* tanh = (exp(2x)-1)/(exp(2x)+1) */
+        inexact = mpfr_mul_2ui (te, x, 1, MPFR_RNDN);  /* 2x */
+        MPFR_ASSERTD(inexact == 0); /* see FIXME above */
+        /* since x > 0, we can only have an overflow */
+        mpfr_exp (te, te, MPFR_RNDN);        /* exp(2x) */
+        if (MPFR_UNLIKELY (MPFR_IS_INF (te)))
           {
-            inexact = -inexact;
-            mpfr_nexttozero (y);
+          set_one:
+            inexact = MPFR_FROM_SIGN_TO_INT (sign);
+            mpfr_set4 (y, __gmpfr_one, MPFR_RNDN, sign);
+            if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG_SIGN (sign)))
+              {
+                inexact = -inexact;
+                mpfr_nexttozero (y);
+              }
+            break;
+          }
+        d = MPFR_GET_EXP (te);               /* For Error calculation */
+        mpfr_add_ui (t, te, 1, MPFR_RNDD);   /* exp(2x) + 1 */
+        mpfr_sub_ui (te, te, 1, MPFR_RNDU);  /* exp(2x) - 1 */
+        d = d - MPFR_GET_EXP (te);
+        mpfr_div (t, te, t, MPFR_RNDN);      /* (exp(2x)-1)/(exp(2x)+1) */
+
+        /* Calculation of the error, see algorithms.tex; the current value
+           of d is k in algorithms.tex. */
+        d = MAX(3, d + 1);  /* d = exponent in 2^(max(3,k+1)) */
+        err = Nt - (d + 1);
+
+        /* The inequality is the condition max(3,k+1) <= floor(p/2). */
+        if (MPFR_LIKELY (d <= Nt / 2 &&
+                         MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
+          {
+            inexact = mpfr_set4 (y, t, rnd_mode, sign);
+            break;
           }
-        break;
-      }
-      d = MPFR_GET_EXP (te);              /* For Error calculation */
-      mpfr_add_ui (t, te, 1, MPFR_RNDD);   /* exp(2x) + 1*/
-      mpfr_sub_ui (te, te, 1, MPFR_RNDU);  /* exp(2x) - 1*/
-      d = d - MPFR_GET_EXP (te);
-      mpfr_div (t, te, t, MPFR_RNDN);      /* (exp(2x)-1)/(exp(2x)+1)*/
-
-      /* Calculation of the error */
-      d = MAX(3, d + 1);
-      err = Nt - (d + 1);
-
-      if (MPFR_LIKELY ((d <= Nt / 2) && MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
-        {
-          inexact = mpfr_set4 (y, t, rnd_mode, sign);
-          break;
-        }
 
-      /* if t=1, we still can round since |sinh(x)| < 1 */
-      if (MPFR_GET_EXP (t) == 1)
-        goto set_one;
+        /* if t=1, we still can round since |sinh(x)| < 1 */
+        if (MPFR_GET_EXP (t) == 1)
+          goto set_one;
 
-      /* Actualisation of the precision */
-      MPFR_ZIV_NEXT (loop, Nt);
-      MPFR_GROUP_REPREC_2 (group, Nt, t, te);
-    }
+        /* Actualisation of the precision */
+        MPFR_ZIV_NEXT (loop, Nt);
+        MPFR_GROUP_REPREC_2 (group, Nt, t, te);
+      }
     MPFR_ZIV_FREE (loop);
     MPFR_GROUP_CLEAR (group);
   }
diff --git a/tests/tdiv.c b/tests/tdiv.c
index 077347e66..0b16cd3a4 100644
--- a/tests/tdiv.c
+++ b/tests/tdiv.c
@@ -1457,6 +1457,30 @@ test_20170105 (void)
   mpfr_clears (x, y, z, t, (mpfr_ptr) 0);
 }
 
+/* The real cause of the mpfr_ttanh failure from the non-regression test
+   added in tests/ttanh.c@11993 was actually due to a bug in mpfr_div, as
+   this can be seen by comparing the logs of the 3.1 branch and the trunk
+   r11992 with MPFR_LOG_ALL=1 MPFR_LOG_PREC=50 on this particular test
+   (this was noticed because adding this test to the 3.1 branch did not
+   yield a failure like in the trunk, though the mpfr_ttanh code did not
+   change until r11993). This was the bug actually fixed in r12002.
+*/
+static void
+test_20171219 (void)
+{
+  mpfr_t x, y, z, t;
+
+  mpfr_inits2 (126, x, y, z, t, (mpfr_ptr) 0);
+  mpfr_set_str_binary (x, "0.111000000000000111100000000000011110000000000001111000000000000111100000000000011110000000000001111000000000000111100000000000E1");
+  /* x = 36347266450315671364380109803814927 / 2^114 */
+  mpfr_add_ui (y, x, 2, MPFR_RNDN);
+  /* y = 77885641318594292392624080437575695 / 2^114 */
+  mpfr_div (z, x, y, MPFR_RNDN);
+  mpfr_set_ui_2exp (t, 3823, -13, MPFR_RNDN);
+  MPFR_ASSERTN (mpfr_equal_p (z, t));
+  mpfr_clears (x, y, z, t, (mpfr_ptr) 0);
+}
+
 #if !defined(MPFR_GENERIC_ABI) && GMP_NUMB_BITS == 64
 /* exercise mpfr_div2_approx */
 static void
@@ -1479,11 +1503,28 @@ test_mpfr_div2_approx (unsigned long n)
 }
 #endif
 
+/* bug found in ttan with GMP_CHECK_RANDOMIZE=1514257254 */
+static void
+bug20171218 (void)
+{
+  mpfr_t s, c;
+  mpfr_init2 (s, 124);
+  mpfr_init2 (c, 124);
+  mpfr_set_str_binary (s, "-0.1110000111100001111000011110000111100001111000011110000111100001111000011110000111100001111000011110000111100001111000011110E0");
+  mpfr_set_str_binary (c, "0.1111000011110000111100001111000011110000111100001111000011110000111100001111000011110000111100001111000011110000111100001111E-1");
+  mpfr_div (c, s, c, MPFR_RNDN);
+  mpfr_set_str_binary (s, "-1.111000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000");
+  MPFR_ASSERTN(mpfr_equal_p (c, s));
+  mpfr_clear (s);
+  mpfr_clear (c);
+}
+
 int
 main (int argc, char *argv[])
 {
   tests_start_mpfr ();
 
+  bug20171218 ();
   testall_rndf (9);
   test_20170105 ();
   check_inexact ();
@@ -1509,6 +1550,7 @@ main (int argc, char *argv[])
   test_20151023 ();
   test_20160831 ();
   test_20170104 ();
+  test_20171219 ();
   test_generic (MPFR_PREC_MIN, 800, 50);
   test_bad ();
   test_extreme ();
diff --git a/tests/texp2.c b/tests/texp2.c
index d568d88e0..e22ed58c4 100644
--- a/tests/texp2.c
+++ b/tests/texp2.c
@@ -215,6 +215,70 @@ overflowed_exp2_0 (void)
   mpfr_clear (y);
 }
 
+static void
+bug20171218 (void)
+{
+  mpfr_t x, y, z;
+  mpfr_exp_t emin;
+  int inex, i;
+  mpfr_flags_t flags, ex_flags;
+
+  emin = mpfr_get_emin ();
+
+  mpfr_init2 (x, 228);
+  mpfr_init2 (y, 11);
+  mpfr_init2 (z, 11);
+  mpfr_set_str_binary (x, "-0.110111010100001100000000000000111001100101011011101110101011000011011011001101111111110100110001110100111000111101010010100010001101100001010111101110100010000101011111001101011000011101000000001010001011110011110101010111000000E17");
+  mpfr_set_emin (-113285);
+  mpfr_clear_flags ();
+  inex = mpfr_exp2 (y, x, MPFR_RNDA);
+  /* exact result is 0.11111111111110110000001011...E-113286, which rounded away
+     gives 0.10000000000E-113285, i.e., no underflow (after rounding) */
+  mpfr_set_str_binary (z, "0.10000000000E-113285");
+  MPFR_ASSERTN(mpfr_equal_p (y, z));
+  MPFR_ASSERTN(inex > 0);
+  MPFR_ASSERTN(mpfr_inexflag_p ());
+  MPFR_ASSERTN(!mpfr_underflow_p ());
+  mpfr_clear (x);
+  mpfr_clear (y);
+  mpfr_clear (z);
+
+  mpfr_init2 (x, 256);
+  mpfr_init2 (y, 8);
+  mpfr_init2 (z, 8);
+
+  for (i = 0; i < 3; i++)
+    {
+      mpfr_set_emin (i == 0 ? -17 : i == 1 ? emin : MPFR_EMIN_MIN);
+      mpfr_set_exp_t (x, __gmpfr_emin - 2, MPFR_RNDN);
+      mpfr_nextabove (x);
+      mpfr_set_ui_2exp (z, 1, __gmpfr_emin - 1, MPFR_RNDN);
+      ex_flags = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
+      mpfr_clear_flags ();
+      inex = mpfr_exp2 (y, x, MPFR_RNDN);
+      flags = __gmpfr_flags;
+      if (! (flags == ex_flags && SAME_SIGN (inex, 1) && mpfr_equal_p (y, z)))
+        {
+          printf ("Error in bug20171218 for i=%d\n", i);
+          printf ("Expected ");
+          mpfr_dump (z);
+          printf ("  with inex = 1 and flags:");
+          flags_out (ex_flags);
+          printf ("Got      ");
+          mpfr_dump (y);
+          printf ("  with inex = %d and flags:", inex);
+          flags_out (flags);
+          exit (1);
+        }
+    }
+
+  mpfr_clear (x);
+  mpfr_clear (y);
+  mpfr_clear (z);
+
+  mpfr_set_emin (emin);
+}
+
 int
 main (int argc, char *argv[])
 {
@@ -223,6 +287,7 @@ main (int argc, char *argv[])
 
   tests_start_mpfr ();
 
+  bug20171218 ();
   special_overflow ();
   emax_m_eps ();
   exp_range ();
diff --git a/tests/tfma.c b/tests/tfma.c
index aa5c894a3..357781665 100644
--- a/tests/tfma.c
+++ b/tests/tfma.c
@@ -485,6 +485,49 @@ bug20101018 (void)
   mpfr_clear (u);
 }
 
+/* bug found with GMP_CHECK_RANDOMIZE=1514407719 */
+static void
+bug20171219 (void)
+{
+  mpfr_t x, y, z, t, u;
+  int i;
+
+  mpfr_init2 (x, 60);
+  mpfr_init2 (y, 60);
+  mpfr_init2 (z, 60);
+  mpfr_init2 (t, 68);
+  mpfr_init2 (u, 68);
+
+  mpfr_set_ui (x, 1, MPFR_RNDN);
+  mpfr_set_ui (y, 1, MPFR_RNDN);
+  mpfr_set_ui (z, 1, MPFR_RNDN);
+  mpfr_set_ui (t, 2, MPFR_RNDN);
+  i = mpfr_fma (u, x, y, z, MPFR_RNDA);
+  if (! mpfr_equal_p (u, t))
+    {
+      printf ("Wrong result in bug20171219 (a)\n");
+      printf ("Expected ");
+      mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN);
+      printf ("\nGot      ");
+      mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN);
+      printf ("\n");
+      exit (1);
+    }
+  if (i != 0)
+    {
+      printf ("Wrong ternary value in bug20171219\n");
+      printf ("Expected 0\n");
+      printf ("Got      %d\n", i);
+      exit (1);
+    }
+
+  mpfr_clear (x);
+  mpfr_clear (y);
+  mpfr_clear (z);
+  mpfr_clear (t);
+  mpfr_clear (u);
+}
+
 int
 main (int argc, char *argv[])
 {
@@ -496,6 +539,7 @@ main (int argc, char *argv[])
   emin = mpfr_get_emin ();
   emax = mpfr_get_emax ();
 
+  bug20171219 ();
   bug20101018 ();
 
   mpfr_init (x);
diff --git a/tests/tmul_d.c b/tests/tmul_d.c
index d29976146..42787b014 100644
--- a/tests/tmul_d.c
+++ b/tests/tmul_d.c
@@ -70,6 +70,24 @@ check_nans (void)
   mpfr_clear (y);
 }
 
+static void
+bug20171218 (void)
+{
+  mpfr_t a, b;
+  mpfr_init2 (a, 22);
+  mpfr_init2 (b, 1312);
+  mpfr_set_str (b, "1.ffffffffffffffffffc3e0007ffe000700fff8001fff800000000007e0fffe0007fffffff00001f800000003ffffffe1fc00000003fffe7c03ffffffffffffffe000000fffffffff83f0000007ffc03ffffffffc0007ff0000ffcfffffe00000f80000000003c007ffffffff3ff807ffffffffff000000000000001fc000fffe000600000ff0003ffe00fffffffc00001ffc0fffffffffff00000807fe03ffffffc01ffe", 16, MPFR_RNDN);
+  mpfr_mul_d (a, b, 0.5, MPFR_RNDF);
+  /* a should be 1 or nextbelow(1) */
+  if (mpfr_cmp_ui (a, 1) != 0)
+    {
+      mpfr_nextabove (a);
+      MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0);
+    }
+  mpfr_clear (a);
+  mpfr_clear (b);
+}
+
 #define TEST_FUNCTION mpfr_mul_d
 #define DOUBLE_ARG2
 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 1, RANDS)
@@ -84,6 +102,8 @@ main (void)
 
   tests_start_mpfr ();
 
+  bug20171218 ();
+
   /* check with enough precision */
   mpfr_init2 (x, IEEE_DBL_MANT_DIG);
   mpfr_init2 (y, IEEE_DBL_MANT_DIG);
diff --git a/tests/ttan.c b/tests/ttan.c
index 557168231..393acfe60 100644
--- a/tests/ttan.c
+++ b/tests/ttan.c
@@ -92,6 +92,25 @@ check_nans (void)
   mpfr_clear (y);
 }
 
+/* This test was generated from bad_cases, with GMP_CHECK_RANDOMIZE=1514257254.
+   For the x value below, we have tan(x) = -1.875 + epsilon,
+   thus with RNDZ it should be rounded to -1.c. */
+static void
+bug20171218 (void)
+{
+  mpfr_t x, y, z;
+  mpfr_init2 (x, 804);
+  mpfr_init2 (y, 4);
+  mpfr_init2 (z, 4);
+  mpfr_set_str (x, "-1.14b1dd5f90ce0eded2a8f59c05e72daf7cc4c78f5075d73246fa420e2c026291d9377e67e7f54e925c4aed39c7b2f917424033c8612f00a19821890558d6c2cef9c60f4ad2c3b061ed53a1709dc1ec8e139627c2119c36d7ebdff0d715e559b47f740c534", 16, MPFR_RNDN);
+  mpfr_tan (y, x, MPFR_RNDZ);
+  mpfr_set_str (z, "-1.c", 16, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_equal_p (y, z));
+  mpfr_clear (x);
+  mpfr_clear (y);
+  mpfr_clear (z);
+}
+
 int
 main (int argc, char *argv[])
 {
@@ -102,6 +121,7 @@ main (int argc, char *argv[])
 
   tests_start_mpfr ();
 
+  bug20171218 ();
   check_nans ();
 
   mpfr_init (x);
diff --git a/tests/ttanh.c b/tests/ttanh.c
index 90798d3ff..e2628a2e3 100644
--- a/tests/ttanh.c
+++ b/tests/ttanh.c
@@ -116,11 +116,31 @@ special_overflow (void)
   mpfr_clear (x);
 }
 
+/* This test was generated from bad_cases, with input y=-7.778@-1 = -3823/8192.
+   For the x value below, we have atanh(y) < x, thus since tanh() is increasing,
+   y < tanh(x), and thus tanh(x) rounded towards zero should give -3822/8192. */
+static void
+bug20171218 (void)
+{
+  mpfr_t x, y, z;
+  mpfr_init2 (x, 813);
+  mpfr_init2 (y, 12);
+  mpfr_init2 (z, 12);
+  mpfr_set_str (x, "-8.17cd20bfc17ae00935dc3abad8e17ab43d3ef7740c320798eefb93191f4a62dba9a2daa5efb6eace21130abd87e3ee2eadd2ad8ddae883d2f2db5dee1ac7ce3c59d16eca09e2ca3f21dc2a0386c037a0d3972e62d5b6e82446032020705553c566b1df24f40@-1", 16, MPFR_RNDN);
+  mpfr_tanh (y, x, MPFR_RNDZ);
+  mpfr_set_str (z, "-7.770@-1", 16, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_equal_p (y, z));
+  mpfr_clear (x);
+  mpfr_clear (y);
+  mpfr_clear (z);
+}
+
 int
 main (int argc, char *argv[])
 {
   tests_start_mpfr ();
 
+  bug20171218 ();
   special_overflow ();
   special ();