Submitted By:            Bruce Dubbs <bdubbs_at_linuxfromscratch_dot_org>
Date:                    2016-02-16
Initial Package Version: 3.1.3
Upstream Status:         Already in upstream patch repo
Origin:                  Upstream http://www.mpfr.org/mpfr-current/allpatches
                         as of 2015-02-15.
Description:             This patch contains bug fixes identified upstream.

--- a/doc/mpfr.info	2015-06-19 21:55:53.000000000 +0200
+++ b/doc/mpfr.info	2016-02-16 14:55:39.806391028 +0100
@@ -1,4 +1,4 @@
-This is mpfr.info, produced by makeinfo version 5.2 from mpfr.texi.
+This is mpfr.info, produced by makeinfo version 6.0 from mpfr.texi.
 
 This manual documents how to install and use the Multiple Precision
 Floating-Point Reliable Library, version 3.1.3.
@@ -55,7 +55,7 @@
 MPFR Copying Conditions
 ***********************
 
-The GNU MPFR library (or MPFR for short) is "free"; this means that
+The GNU MPFR library (or MPFR for short) is “free”; this means that
 everyone is free to use it and free to redistribute it on a free basis.
 The library is not in the public domain; it is copyrighted and there are
 restrictions on its distribution, but these restrictions are designed to
@@ -418,7 +418,7 @@
 4.2 Nomenclature and Types
 ==========================
 
-A "floating-point number", or "float" for short, is an arbitrary
+A “floating-point number”, or “float” for short, is an arbitrary
 precision significand (also called mantissa) with a limited precision
 exponent.  The C data type for such objects is ‘mpfr_t’ (internally
 defined as a one-element array of a structure, and ‘mpfr_ptr’ is the C
@@ -432,7 +432,7 @@
 to the other functions supported by MPFR. Unless documented otherwise,
 the sign bit of a NaN is unspecified.
 
-The "precision" is the number of bits used to represent the significand
+The “precision” is the number of bits used to represent the significand
 of a floating-point number; the corresponding C data type is
 ‘mpfr_prec_t’.  The precision can be any integer between ‘MPFR_PREC_MIN’
 and ‘MPFR_PREC_MAX’.  In the current implementation, ‘MPFR_PREC_MIN’ is
@@ -446,7 +446,7 @@
 may abort, crash or have undefined behavior (depending on your C
 implementation).
 
-The "rounding mode" specifies the way to round the result of a
+The “rounding mode” specifies the way to round the result of a
 floating-point operation, in case the exact result can not be
 represented exactly in the destination significand; the corresponding C
 data type is ‘mpfr_rnd_t’.
@@ -499,14 +499,14 @@
 representable numbers, it is rounded to the one with the least
 significant bit set to zero.  For example, the number 2.5, which is
 represented by (10.1) in binary, is rounded to (10.0)=2 with a precision
-of two bits, and not to (11.0)=3.  This rule avoids the "drift"
+of two bits, and not to (11.0)=3.  This rule avoids the “drift”
 phenomenon mentioned by Knuth in volume 2 of The Art of Computer
 Programming (Section 4.2.2).
 
    Most MPFR functions take as first argument the destination variable,
 as second and following arguments the input variables, as last argument
 a rounding mode, and have a return value of type ‘int’, called the
-"ternary value".  The value stored in the destination variable is
+“ternary value”.  The value stored in the destination variable is
 correctly rounded, i.e., MPFR behaves as if it computed the result with
 an infinite precision, then rounded it to the precision of this
 variable.  The input variables are regarded as exact (in particular,
@@ -572,15 +572,18 @@
    When the input point is in the closure of the domain of the
 mathematical function and an input argument is +0 (resp. −0), one
 considers the limit when the corresponding argument approaches 0 from
-above (resp. below).  If the limit is not defined (e.g., ‘mpfr_log’ on
-−0), the behavior is specified in the description of the MPFR function.
+above (resp. below), if possible.  If the limit is not defined (e.g.,
+‘mpfr_sqrt’ and ‘mpfr_log’ on −0), the behavior is specified in the
+description of the MPFR function, but must be consistent with the rule
+from the above paragraph (e.g., ‘mpfr_log’ on ±0 gives −Inf).
 
    When the result is equal to 0, its sign is determined by considering
 the limit as if the input point were not in the domain: If one
 approaches 0 from above (resp. below), the result is +0 (resp. −0); for
-example, ‘mpfr_sin’ on +0 gives +0.  In the other cases, the sign is
-specified in the description of the MPFR function; for example
-‘mpfr_max’ on −0 and +0 gives +0.
+example, ‘mpfr_sin’ on −0 gives −0 and ‘mpfr_acos’ on 1 gives +0 (in all
+rounding modes).  In the other cases, the sign is specified in the
+description of the MPFR function; for example ‘mpfr_max’ on −0 and +0
+gives +0.
 
    When the input point is not in the closure of the domain of the
 function, the result is NaN. Example: ‘mpfr_sqrt’ on −17 gives NaN.
@@ -590,8 +593,8 @@
 numbers; such a case is always explicitly specified in *note MPFR
 Interface::.  Example: ‘mpfr_hypot’ on (NaN,0) gives NaN, but
 ‘mpfr_hypot’ on (NaN,+Inf) gives +Inf (as specified in *note Special
-Functions::), since for any finite input X, ‘mpfr_hypot’ on (X,+Inf)
-gives +Inf.
+Functions::), since for any finite or infinite input X, ‘mpfr_hypot’ on
+(X,+Inf) gives +Inf.
 
 
 File: mpfr.info,  Node: Exceptions,  Next: Memory Handling,  Prev: Floating-Point Values on Special Numbers,  Up: MPFR Basics
@@ -1253,8 +1256,9 @@
           mpfr_rnd_t RND)
  -- Function: int mpfr_add_q (mpfr_t ROP, mpfr_t OP1, mpq_t OP2,
           mpfr_rnd_t RND)
-     Set ROP to OP1 + OP2 rounded in the direction RND.  For types
-     having no signed zero, it is considered unsigned (i.e., (+0) + 0 =
+     Set ROP to OP1 + OP2 rounded in the direction RND.  The IEEE-754
+     rules are used, in particular for signed zeros.  But for types
+     having no signed zeros, 0 is considered unsigned (i.e., (+0) + 0 =
      (+0) and (−0) + 0 = (−0)).  The ‘mpfr_add_d’ function assumes that
      the radix of the ‘double’ type is a power of 2, with a precision at
      most that declared by the C implementation (macro
@@ -1280,8 +1284,9 @@
           mpfr_rnd_t RND)
  -- Function: int mpfr_sub_q (mpfr_t ROP, mpfr_t OP1, mpq_t OP2,
           mpfr_rnd_t RND)
-     Set ROP to OP1 - OP2 rounded in the direction RND.  For types
-     having no signed zero, it is considered unsigned (i.e., (+0) − 0 =
+     Set ROP to OP1 - OP2 rounded in the direction RND.  The IEEE-754
+     rules are used, in particular for signed zeros.  But for types
+     having no signed zeros, 0 is considered unsigned (i.e., (+0) − 0 =
      (+0), (−0) − 0 = (−0), 0 − (+0) = (−0) and 0 − (−0) = (+0)).  The
      same restrictions than for ‘mpfr_add_d’ apply to ‘mpfr_d_sub’ and
      ‘mpfr_sub_d’.
@@ -1300,7 +1305,7 @@
           mpfr_rnd_t RND)
      Set ROP to OP1 times OP2 rounded in the direction RND.  When a
      result is zero, its sign is the product of the signs of the
-     operands (for types having no signed zero, it is considered
+     operands (for types having no signed zeros, 0 is considered
      positive).  The same restrictions than for ‘mpfr_add_d’ apply to
      ‘mpfr_mul_d’.
 
@@ -1327,21 +1332,24 @@
           mpfr_rnd_t RND)
      Set ROP to OP1/OP2 rounded in the direction RND.  When a result is
      zero, its sign is the product of the signs of the operands (for
-     types having no signed zero, it is considered positive).  The same
+     types having no signed zeros, 0 is considered positive).  The same
      restrictions than for ‘mpfr_add_d’ apply to ‘mpfr_d_div’ and
      ‘mpfr_div_d’.
 
  -- Function: int mpfr_sqrt (mpfr_t ROP, mpfr_t OP, mpfr_rnd_t RND)
  -- Function: int mpfr_sqrt_ui (mpfr_t ROP, unsigned long int OP,
           mpfr_rnd_t RND)
-     Set ROP to the square root of OP rounded in the direction RND (set
-     ROP to −0 if OP is −0, to be consistent with the IEEE 754
-     standard).  Set ROP to NaN if OP is negative.
+     Set ROP to the square root of OP rounded in the direction RND.  Set
+     ROP to −0 if OP is −0, to be consistent with the IEEE 754 standard.
+     Set ROP to NaN if OP is negative.
 
  -- Function: int mpfr_rec_sqrt (mpfr_t ROP, mpfr_t OP, mpfr_rnd_t RND)
      Set ROP to the reciprocal square root of OP rounded in the
      direction RND.  Set ROP to +Inf if OP is ±0, +0 if OP is +Inf, and
-     NaN if OP is negative.
+     NaN if OP is negative.  Warning!  Therefore the result on −0 is
+     different from the one of the rSqrt function recommended by the
+     IEEE 754-2008 standard (Section 9.2.1), which is −Inf instead of
+     +Inf.
 
  -- Function: int mpfr_cbrt (mpfr_t ROP, mpfr_t OP, mpfr_rnd_t RND)
  -- Function: int mpfr_root (mpfr_t ROP, mpfr_t OP, unsigned long int K,
@@ -1515,8 +1523,10 @@
  -- Function: int mpfr_log2 (mpfr_t ROP, mpfr_t OP, mpfr_rnd_t RND)
  -- Function: int mpfr_log10 (mpfr_t ROP, mpfr_t OP, mpfr_rnd_t RND)
      Set ROP to the natural logarithm of OP, log2(OP) or log10(OP),
-     respectively, rounded in the direction RND.  Set ROP to −Inf if OP
-     is −0 (i.e., the sign of the zero has no influence on the result).
+     respectively, rounded in the direction RND.  Set ROP to +0 if OP is
+     1 (in all rounding modes), for consistency with the ISO C99 and
+     IEEE 754-2008 standards.  Set ROP to −Inf if OP is ±0 (i.e., the
+     sign of the zero has no influence on the result).
 
  -- Function: int mpfr_exp (mpfr_t ROP, mpfr_t OP, mpfr_rnd_t RND)
  -- Function: int mpfr_exp2 (mpfr_t ROP, mpfr_t OP, mpfr_rnd_t RND)
@@ -1649,17 +1659,21 @@
 
  -- Function: int mpfr_lngamma (mpfr_t ROP, mpfr_t OP, mpfr_rnd_t RND)
      Set ROP to the value of the logarithm of the Gamma function on OP,
-     rounded in the direction RND.  When −2K−1 <= OP <= −2K, K being a
-     non-negative integer, ROP is set to NaN. See also ‘mpfr_lgamma’.
+     rounded in the direction RND.  When OP is 1 or 2, set ROP to +0 (in
+     all rounding modes).  When OP is an infinity or a nonpositive
+     integer, set ROP to +Inf, following the general rules on special
+     values.  When −2K−1 < OP < −2K, K being a nonnegative integer, set
+     ROP to NaN.  See also ‘mpfr_lgamma’.
 
  -- Function: int mpfr_lgamma (mpfr_t ROP, int *SIGNP, mpfr_t OP,
           mpfr_rnd_t RND)
      Set ROP to the value of the logarithm of the absolute value of the
      Gamma function on OP, rounded in the direction RND.  The sign (1 or
      −1) of Gamma(OP) is returned in the object pointed to by SIGNP.
-     When OP is an infinity or a non-positive integer, set ROP to +Inf.
-     When OP is NaN, −Inf or a negative integer, *SIGNP is undefined,
-     and when OP is ±0, *SIGNP is the sign of the zero.
+     When OP is 1 or 2, set ROP to +0 (in all rounding modes).  When OP
+     is an infinity or a nonpositive integer, set ROP to +Inf.  When OP
+     is NaN, −Inf or a negative integer, *SIGNP is undefined, and when
+     OP is ±0, *SIGNP is the sign of the zero.
 
  -- Function: int mpfr_digamma (mpfr_t ROP, mpfr_t OP, mpfr_rnd_t RND)
      Set ROP to the value of the Digamma (sometimes also called Psi)
@@ -1703,7 +1717,10 @@
  -- Function: int mpfr_fms (mpfr_t ROP, mpfr_t OP1, mpfr_t OP2, mpfr_t
           OP3, mpfr_rnd_t RND)
      Set ROP to (OP1 times OP2) + OP3 (resp. (OP1 times OP2) - OP3)
-     rounded in the direction RND.
+     rounded in the direction RND.  Concerning special values (signed
+     zeros, infinities, NaN), these functions behave like a
+     multiplication followed by a separate addition or subtraction.
+     That is, the fused operation matters only for rounding.
 
  -- Function: int mpfr_agm (mpfr_t ROP, mpfr_t OP1, mpfr_t OP2,
           mpfr_rnd_t RND)
@@ -1717,9 +1734,10 @@
           RND)
      Set ROP to the Euclidean norm of X and Y, i.e., the square root of
      the sum of the squares of X and Y, rounded in the direction RND.
-     Special values are handled as described in Section F.9.4.3 of the
-     ISO C99 and IEEE 754-2008 standards: If X or Y is an infinity, then
-     +Inf is returned in ROP, even if the other number is NaN.
+     Special values are handled as described in the ISO C99 (Section
+     F.9.4.3) and IEEE 754-2008 (Section 9.2.1) standards: If X or Y is
+     an infinity, then +Inf is returned in ROP, even if the other number
+     is NaN.
 
  -- Function: int mpfr_ai (mpfr_t ROP, mpfr_t X, mpfr_rnd_t RND)
      Set ROP to the value of the Airy function Ai on X, rounded in the
@@ -2670,7 +2688,7 @@
 5.16 Internals
 ==============
 
-A "limb" means the part of a multi-precision number that fits in a
+A “limb” means the part of a multi-precision number that fits in a
 single word.  Usually a limb contains 32 or 64 bits.  The C data type
 for a limb is ‘mp_limb_t’.
 
@@ -3140,7 +3158,7 @@
   0. PREAMBLE
 
      The purpose of this License is to make a manual, textbook, or other
-     functional and useful document "free" in the sense of freedom: to
+     functional and useful document “free” in the sense of freedom: to
      assure everyone the effective freedom to copy and redistribute it,
      with or without modifying it, either commercially or
      noncommercially.  Secondarily, this License preserves for the
@@ -3655,9 +3673,9 @@
 * Menu:
 
 * mpfr_abs:                              Basic Arithmetic Functions.
-                                                              (line 160)
-* mpfr_acos:                             Special Functions.   (line  51)
-* mpfr_acosh:                            Special Functions.   (line 115)
+                                                              (line 165)
+* mpfr_acos:                             Special Functions.   (line  53)
+* mpfr_acosh:                            Special Functions.   (line 117)
 * mpfr_add:                              Basic Arithmetic Functions.
                                                               (line   6)
 * mpfr_add_d:                            Basic Arithmetic Functions.
@@ -3670,15 +3688,15 @@
                                                               (line   8)
 * mpfr_add_z:                            Basic Arithmetic Functions.
                                                               (line  14)
-* mpfr_agm:                              Special Functions.   (line 210)
-* mpfr_ai:                               Special Functions.   (line 226)
-* mpfr_asin:                             Special Functions.   (line  52)
-* mpfr_asinh:                            Special Functions.   (line 116)
+* mpfr_agm:                              Special Functions.   (line 219)
+* mpfr_ai:                               Special Functions.   (line 236)
+* mpfr_asin:                             Special Functions.   (line  54)
+* mpfr_asinh:                            Special Functions.   (line 118)
 * mpfr_asprintf:                         Formatted Output Functions.
                                                               (line 193)
-* mpfr_atan:                             Special Functions.   (line  53)
-* mpfr_atan2:                            Special Functions.   (line  63)
-* mpfr_atanh:                            Special Functions.   (line 117)
+* mpfr_atan:                             Special Functions.   (line  55)
+* mpfr_atan2:                            Special Functions.   (line  65)
+* mpfr_atanh:                            Special Functions.   (line 119)
 * mpfr_buildopt_decimal_p:               Miscellaneous Functions.
                                                               (line 162)
 * mpfr_buildopt_gmpinternals_p:          Miscellaneous Functions.
@@ -3690,7 +3708,7 @@
 * mpfr_can_round:                        Rounding Related Functions.
                                                               (line  39)
 * mpfr_cbrt:                             Basic Arithmetic Functions.
-                                                              (line 108)
+                                                              (line 113)
 * mpfr_ceil:                             Integer Related Functions.
                                                               (line   7)
 * mpfr_check_range:                      Exception Related Functions.
@@ -3735,18 +3753,18 @@
                                                               (line  27)
 * mpfr_cmp_z:                            Comparison Functions.
                                                               (line  11)
-* mpfr_const_catalan:                    Special Functions.   (line 237)
-* mpfr_const_euler:                      Special Functions.   (line 236)
-* mpfr_const_log2:                       Special Functions.   (line 234)
-* mpfr_const_pi:                         Special Functions.   (line 235)
+* mpfr_const_catalan:                    Special Functions.   (line 247)
+* mpfr_const_euler:                      Special Functions.   (line 246)
+* mpfr_const_log2:                       Special Functions.   (line 244)
+* mpfr_const_pi:                         Special Functions.   (line 245)
 * mpfr_copysign:                         Miscellaneous Functions.
                                                               (line 109)
-* mpfr_cos:                              Special Functions.   (line  29)
-* mpfr_cosh:                             Special Functions.   (line  95)
-* mpfr_cot:                              Special Functions.   (line  47)
-* mpfr_coth:                             Special Functions.   (line 111)
-* mpfr_csc:                              Special Functions.   (line  46)
-* mpfr_csch:                             Special Functions.   (line 110)
+* mpfr_cos:                              Special Functions.   (line  31)
+* mpfr_cosh:                             Special Functions.   (line  97)
+* mpfr_cot:                              Special Functions.   (line  49)
+* mpfr_coth:                             Special Functions.   (line 113)
+* mpfr_csc:                              Special Functions.   (line  48)
+* mpfr_csch:                             Special Functions.   (line 112)
 * mpfr_custom_get_exp:                   Custom Interface.    (line  75)
 * mpfr_custom_get_kind:                  Custom Interface.    (line  65)
 * mpfr_custom_get_significand:           Custom Interface.    (line  70)
@@ -3756,47 +3774,47 @@
 * mpfr_custom_move:                      Custom Interface.    (line  82)
 * MPFR_DECL_INIT:                        Initialization Functions.
                                                               (line  74)
-* mpfr_digamma:                          Special Functions.   (line 166)
+* mpfr_digamma:                          Special Functions.   (line 172)
 * mpfr_dim:                              Basic Arithmetic Functions.
-                                                              (line 166)
+                                                              (line 171)
 * mpfr_div:                              Basic Arithmetic Functions.
-                                                              (line  72)
+                                                              (line  74)
 * mpfr_divby0_p:                         Exception Related Functions.
                                                               (line 134)
 * mpfr_div_2exp:                         Compatibility with MPF.
                                                               (line  49)
 * mpfr_div_2si:                          Basic Arithmetic Functions.
-                                                              (line 181)
+                                                              (line 186)
 * mpfr_div_2ui:                          Basic Arithmetic Functions.
-                                                              (line 179)
+                                                              (line 184)
 * mpfr_div_d:                            Basic Arithmetic Functions.
-                                                              (line  84)
+                                                              (line  86)
 * mpfr_div_q:                            Basic Arithmetic Functions.
-                                                              (line  88)
+                                                              (line  90)
 * mpfr_div_si:                           Basic Arithmetic Functions.
-                                                              (line  80)
+                                                              (line  82)
 * mpfr_div_ui:                           Basic Arithmetic Functions.
-                                                              (line  76)
+                                                              (line  78)
 * mpfr_div_z:                            Basic Arithmetic Functions.
-                                                              (line  86)
+                                                              (line  88)
 * mpfr_d_div:                            Basic Arithmetic Functions.
-                                                              (line  82)
+                                                              (line  84)
 * mpfr_d_sub:                            Basic Arithmetic Functions.
-                                                              (line  35)
-* mpfr_eint:                             Special Functions.   (line 133)
+                                                              (line  36)
+* mpfr_eint:                             Special Functions.   (line 135)
 * mpfr_eq:                               Compatibility with MPF.
                                                               (line  28)
 * mpfr_equal_p:                          Comparison Functions.
                                                               (line  59)
 * mpfr_erangeflag_p:                     Exception Related Functions.
                                                               (line 137)
-* mpfr_erf:                              Special Functions.   (line 177)
-* mpfr_erfc:                             Special Functions.   (line 178)
-* mpfr_exp:                              Special Functions.   (line  23)
-* mpfr_exp10:                            Special Functions.   (line  25)
-* mpfr_exp2:                             Special Functions.   (line  24)
-* mpfr_expm1:                            Special Functions.   (line 129)
-* mpfr_fac_ui:                           Special Functions.   (line 121)
+* mpfr_erf:                              Special Functions.   (line 183)
+* mpfr_erfc:                             Special Functions.   (line 184)
+* mpfr_exp:                              Special Functions.   (line  25)
+* mpfr_exp10:                            Special Functions.   (line  27)
+* mpfr_exp2:                             Special Functions.   (line  26)
+* mpfr_expm1:                            Special Functions.   (line 131)
+* mpfr_fac_ui:                           Special Functions.   (line 123)
 * mpfr_fits_intmax_p:                    Conversion Functions.
                                                               (line 150)
 * mpfr_fits_sint_p:                      Conversion Functions.
@@ -3815,20 +3833,20 @@
                                                               (line 147)
 * mpfr_floor:                            Integer Related Functions.
                                                               (line   8)
-* mpfr_fma:                              Special Functions.   (line 203)
+* mpfr_fma:                              Special Functions.   (line 209)
 * mpfr_fmod:                             Integer Related Functions.
                                                               (line  92)
-* mpfr_fms:                              Special Functions.   (line 205)
+* mpfr_fms:                              Special Functions.   (line 211)
 * mpfr_fprintf:                          Formatted Output Functions.
                                                               (line 157)
 * mpfr_frac:                             Integer Related Functions.
                                                               (line  76)
-* mpfr_free_cache:                       Special Functions.   (line 244)
+* mpfr_free_cache:                       Special Functions.   (line 254)
 * mpfr_free_str:                         Conversion Functions.
                                                               (line 137)
 * mpfr_frexp:                            Conversion Functions.
                                                               (line  45)
-* mpfr_gamma:                            Special Functions.   (line 148)
+* mpfr_gamma:                            Special Functions.   (line 150)
 * mpfr_get_d:                            Conversion Functions.
                                                               (line   7)
 * mpfr_get_decimal64:                    Conversion Functions.
@@ -3887,7 +3905,7 @@
                                                               (line  56)
 * mpfr_greater_p:                        Comparison Functions.
                                                               (line  55)
-* mpfr_hypot:                            Special Functions.   (line 218)
+* mpfr_hypot:                            Special Functions.   (line 227)
 * mpfr_inexflag_p:                       Exception Related Functions.
                                                               (line 136)
 * mpfr_inf_p:                            Comparison Functions.
@@ -3922,21 +3940,21 @@
                                                               (line  31)
 * mpfr_integer_p:                        Integer Related Functions.
                                                               (line 119)
-* mpfr_j0:                               Special Functions.   (line 182)
-* mpfr_j1:                               Special Functions.   (line 183)
-* mpfr_jn:                               Special Functions.   (line 184)
+* mpfr_j0:                               Special Functions.   (line 188)
+* mpfr_j1:                               Special Functions.   (line 189)
+* mpfr_jn:                               Special Functions.   (line 190)
 * mpfr_lessequal_p:                      Comparison Functions.
                                                               (line  58)
 * mpfr_lessgreater_p:                    Comparison Functions.
                                                               (line  64)
 * mpfr_less_p:                           Comparison Functions.
                                                               (line  57)
-* mpfr_lgamma:                           Special Functions.   (line 157)
-* mpfr_li2:                              Special Functions.   (line 143)
-* mpfr_lngamma:                          Special Functions.   (line 152)
+* mpfr_lgamma:                           Special Functions.   (line 162)
+* mpfr_li2:                              Special Functions.   (line 145)
+* mpfr_lngamma:                          Special Functions.   (line 154)
 * mpfr_log:                              Special Functions.   (line  16)
 * mpfr_log10:                            Special Functions.   (line  18)
-* mpfr_log1p:                            Special Functions.   (line 125)
+* mpfr_log1p:                            Special Functions.   (line 127)
 * mpfr_log2:                             Special Functions.   (line  17)
 * mpfr_max:                              Miscellaneous Functions.
                                                               (line  22)
@@ -3947,29 +3965,29 @@
 * mpfr_modf:                             Integer Related Functions.
                                                               (line  82)
 * mpfr_mul:                              Basic Arithmetic Functions.
-                                                              (line  51)
+                                                              (line  53)
 * mpfr_mul_2exp:                         Compatibility with MPF.
                                                               (line  47)
 * mpfr_mul_2si:                          Basic Arithmetic Functions.
-                                                              (line 174)
+                                                              (line 179)
 * mpfr_mul_2ui:                          Basic Arithmetic Functions.
-                                                              (line 172)
+                                                              (line 177)
 * mpfr_mul_d:                            Basic Arithmetic Functions.
-                                                              (line  57)
+                                                              (line  59)
 * mpfr_mul_q:                            Basic Arithmetic Functions.
-                                                              (line  61)
+                                                              (line  63)
 * mpfr_mul_si:                           Basic Arithmetic Functions.
-                                                              (line  55)
+                                                              (line  57)
 * mpfr_mul_ui:                           Basic Arithmetic Functions.
-                                                              (line  53)
+                                                              (line  55)
 * mpfr_mul_z:                            Basic Arithmetic Functions.
-                                                              (line  59)
+                                                              (line  61)
 * mpfr_nanflag_p:                        Exception Related Functions.
                                                               (line 135)
 * mpfr_nan_p:                            Comparison Functions.
                                                               (line  39)
 * mpfr_neg:                              Basic Arithmetic Functions.
-                                                              (line 159)
+                                                              (line 164)
 * mpfr_nextabove:                        Miscellaneous Functions.
                                                               (line  15)
 * mpfr_nextbelow:                        Miscellaneous Functions.
@@ -3983,13 +4001,13 @@
 * mpfr_overflow_p:                       Exception Related Functions.
                                                               (line 133)
 * mpfr_pow:                              Basic Arithmetic Functions.
-                                                              (line 116)
+                                                              (line 121)
 * mpfr_pow_si:                           Basic Arithmetic Functions.
-                                                              (line 120)
+                                                              (line 125)
 * mpfr_pow_ui:                           Basic Arithmetic Functions.
-                                                              (line 118)
+                                                              (line 123)
 * mpfr_pow_z:                            Basic Arithmetic Functions.
-                                                              (line 122)
+                                                              (line 127)
 * mpfr_prec_round:                       Rounding Related Functions.
                                                               (line  13)
 * ‘mpfr_prec_t’:                         Nomenclature and Types.
@@ -3999,7 +4017,7 @@
 * mpfr_print_rnd_mode:                   Rounding Related Functions.
                                                               (line  71)
 * mpfr_rec_sqrt:                         Basic Arithmetic Functions.
-                                                              (line 103)
+                                                              (line 105)
 * mpfr_regular_p:                        Comparison Functions.
                                                               (line  43)
 * mpfr_reldiff:                          Compatibility with MPF.
@@ -4021,11 +4039,11 @@
 * ‘mpfr_rnd_t’:                          Nomenclature and Types.
                                                               (line  34)
 * mpfr_root:                             Basic Arithmetic Functions.
-                                                              (line 109)
+                                                              (line 114)
 * mpfr_round:                            Integer Related Functions.
                                                               (line   9)
-* mpfr_sec:                              Special Functions.   (line  45)
-* mpfr_sech:                             Special Functions.   (line 109)
+* mpfr_sec:                              Special Functions.   (line  47)
+* mpfr_sech:                             Special Functions.   (line 111)
 * mpfr_set:                              Assignment Functions.
                                                               (line   9)
 * mpfr_setsign:                          Miscellaneous Functions.
@@ -4100,57 +4118,57 @@
                                                               (line  49)
 * mpfr_signbit:                          Miscellaneous Functions.
                                                               (line  99)
-* mpfr_sin:                              Special Functions.   (line  30)
-* mpfr_sinh:                             Special Functions.   (line  96)
-* mpfr_sinh_cosh:                        Special Functions.   (line 101)
-* mpfr_sin_cos:                          Special Functions.   (line  35)
+* mpfr_sin:                              Special Functions.   (line  32)
+* mpfr_sinh:                             Special Functions.   (line  98)
+* mpfr_sinh_cosh:                        Special Functions.   (line 103)
+* mpfr_sin_cos:                          Special Functions.   (line  37)
 * mpfr_si_div:                           Basic Arithmetic Functions.
-                                                              (line  78)
+                                                              (line  80)
 * mpfr_si_sub:                           Basic Arithmetic Functions.
-                                                              (line  31)
+                                                              (line  32)
 * mpfr_snprintf:                         Formatted Output Functions.
                                                               (line 180)
 * mpfr_sprintf:                          Formatted Output Functions.
                                                               (line 170)
 * mpfr_sqr:                              Basic Arithmetic Functions.
-                                                              (line  69)
+                                                              (line  71)
 * mpfr_sqrt:                             Basic Arithmetic Functions.
-                                                              (line  96)
+                                                              (line  98)
 * mpfr_sqrt_ui:                          Basic Arithmetic Functions.
-                                                              (line  97)
+                                                              (line  99)
 * mpfr_strtofr:                          Assignment Functions.
                                                               (line  80)
 * mpfr_sub:                              Basic Arithmetic Functions.
-                                                              (line  25)
+                                                              (line  26)
 * mpfr_subnormalize:                     Exception Related Functions.
                                                               (line  60)
 * mpfr_sub_d:                            Basic Arithmetic Functions.
-                                                              (line  37)
+                                                              (line  38)
 * mpfr_sub_q:                            Basic Arithmetic Functions.
-                                                              (line  43)
+                                                              (line  44)
 * mpfr_sub_si:                           Basic Arithmetic Functions.
-                                                              (line  33)
+                                                              (line  34)
 * mpfr_sub_ui:                           Basic Arithmetic Functions.
-                                                              (line  29)
+                                                              (line  30)
 * mpfr_sub_z:                            Basic Arithmetic Functions.
-                                                              (line  41)
-* mpfr_sum:                              Special Functions.   (line 252)
+                                                              (line  42)
+* mpfr_sum:                              Special Functions.   (line 262)
 * mpfr_swap:                             Assignment Functions.
                                                               (line 150)
 * ‘mpfr_t’:                              Nomenclature and Types.
                                                               (line   6)
-* mpfr_tan:                              Special Functions.   (line  31)
-* mpfr_tanh:                             Special Functions.   (line  97)
+* mpfr_tan:                              Special Functions.   (line  33)
+* mpfr_tanh:                             Special Functions.   (line  99)
 * mpfr_trunc:                            Integer Related Functions.
                                                               (line  10)
 * mpfr_ui_div:                           Basic Arithmetic Functions.
-                                                              (line  74)
+                                                              (line  76)
 * mpfr_ui_pow:                           Basic Arithmetic Functions.
-                                                              (line 126)
+                                                              (line 131)
 * mpfr_ui_pow_ui:                        Basic Arithmetic Functions.
-                                                              (line 124)
+                                                              (line 129)
 * mpfr_ui_sub:                           Basic Arithmetic Functions.
-                                                              (line  27)
+                                                              (line  28)
 * mpfr_underflow_p:                      Exception Related Functions.
                                                               (line 132)
 * mpfr_unordered_p:                      Comparison Functions.
@@ -4181,61 +4199,61 @@
                                                               (line 182)
 * mpfr_vsprintf:                         Formatted Output Functions.
                                                               (line 171)
-* mpfr_y0:                               Special Functions.   (line 193)
-* mpfr_y1:                               Special Functions.   (line 194)
-* mpfr_yn:                               Special Functions.   (line 195)
+* mpfr_y0:                               Special Functions.   (line 199)
+* mpfr_y1:                               Special Functions.   (line 200)
+* mpfr_yn:                               Special Functions.   (line 201)
 * mpfr_zero_p:                           Comparison Functions.
                                                               (line  42)
-* mpfr_zeta:                             Special Functions.   (line 171)
-* mpfr_zeta_ui:                          Special Functions.   (line 172)
+* mpfr_zeta:                             Special Functions.   (line 177)
+* mpfr_zeta_ui:                          Special Functions.   (line 178)
 * mpfr_z_sub:                            Basic Arithmetic Functions.
-                                                              (line  39)
+                                                              (line  40)
 
 
 
 Tag Table:
 Node: Top775
 Node: Copying2007
-Node: Introduction to MPFR3766
-Node: Installing MPFR5880
-Node: Reporting Bugs11323
-Node: MPFR Basics13353
-Node: Headers and Libraries13669
-Node: Nomenclature and Types16828
-Node: MPFR Variable Conventions18874
-Node: Rounding Modes20418
-Ref: ternary value21544
-Node: Floating-Point Values on Special Numbers23526
-Node: Exceptions26572
-Node: Memory Handling29749
-Node: MPFR Interface30894
-Node: Initialization Functions33008
-Node: Assignment Functions40318
-Node: Combined Initialization and Assignment Functions49673
-Node: Conversion Functions50974
-Node: Basic Arithmetic Functions60035
-Node: Comparison Functions69200
-Node: Special Functions72687
-Node: Input and Output Functions86672
-Node: Formatted Output Functions88644
-Node: Integer Related Functions98431
-Node: Rounding Related Functions105051
-Node: Miscellaneous Functions108888
-Node: Exception Related Functions117568
-Node: Compatibility with MPF124386
-Node: Custom Interface127127
-Node: Internals131526
-Node: API Compatibility133066
-Node: Type and Macro Changes134995
-Node: Added Functions137844
-Node: Changed Functions141132
-Node: Removed Functions145545
-Node: Other Changes145973
-Node: Contributors147576
-Node: References150219
-Node: GNU Free Documentation License151973
-Node: Concept Index174562
-Node: Function and Type Index180659
+Node: Introduction to MPFR3770
+Node: Installing MPFR5884
+Node: Reporting Bugs11327
+Node: MPFR Basics13357
+Node: Headers and Libraries13673
+Node: Nomenclature and Types16832
+Node: MPFR Variable Conventions18894
+Node: Rounding Modes20438
+Ref: ternary value21568
+Node: Floating-Point Values on Special Numbers23554
+Node: Exceptions26813
+Node: Memory Handling29990
+Node: MPFR Interface31135
+Node: Initialization Functions33249
+Node: Assignment Functions40559
+Node: Combined Initialization and Assignment Functions49914
+Node: Conversion Functions51215
+Node: Basic Arithmetic Functions60276
+Node: Comparison Functions69777
+Node: Special Functions73264
+Node: Input and Output Functions87862
+Node: Formatted Output Functions89834
+Node: Integer Related Functions99621
+Node: Rounding Related Functions106241
+Node: Miscellaneous Functions110078
+Node: Exception Related Functions118758
+Node: Compatibility with MPF125576
+Node: Custom Interface128317
+Node: Internals132716
+Node: API Compatibility134260
+Node: Type and Macro Changes136189
+Node: Added Functions139038
+Node: Changed Functions142326
+Node: Removed Functions146739
+Node: Other Changes147167
+Node: Contributors148770
+Node: References151413
+Node: GNU Free Documentation License153167
+Node: Concept Index175760
+Node: Function and Type Index181857
 
 End Tag Table
 
--- a/doc/mpfr.texi	2015-06-19 21:55:11.000000000 +0200
+++ b/doc/mpfr.texi	2016-02-16 14:55:39.802390950 +0100
@@ -810,13 +810,17 @@
 When the input point is in the closure of the domain of the mathematical
 function and an input argument is +0 (resp.@: @minus{}0), one considers
 the limit when the corresponding argument approaches 0 from above
-(resp.@: below). If the limit is not defined (e.g., @code{mpfr_log} on
-@minus{}0), the behavior is specified in the description of the MPFR function.
+(resp.@: below), if possible. If the limit is not defined (e.g.,
+@code{mpfr_sqrt} and @code{mpfr_log} on @minus{}0), the behavior is
+specified in the description of the MPFR function, but must be consistent
+with the rule from the above paragraph (e.g., @code{mpfr_log} on @pom{}0
+gives @minus{}Inf).
 
 When the result is equal to 0, its sign is determined by considering the
 limit as if the input point were not in the domain: If one approaches 0
 from above (resp.@: below), the result is +0 (resp.@: @minus{}0);
-for example, @code{mpfr_sin} on +0 gives +0.
+for example, @code{mpfr_sin} on @minus{}0 gives @minus{}0 and
+@code{mpfr_acos} on 1 gives +0 (in all rounding modes).
 In the other cases, the sign is specified in the description of the MPFR
 function; for example @code{mpfr_max} on @minus{}0 and +0 gives +0.
 
@@ -832,8 +836,8 @@
 @c that advantages in practice), like for any bug fix.
 Example: @code{mpfr_hypot} on (NaN,0) gives NaN, but @code{mpfr_hypot}
 on (NaN,+Inf) gives +Inf (as specified in @ref{Special Functions}),
-since for any finite input @var{x}, @code{mpfr_hypot} on (@var{x},+Inf)
-gives +Inf.
+since for any finite or infinite input @var{x}, @code{mpfr_hypot} on
+(@var{x},+Inf) gives +Inf.
 
 @node Exceptions, Memory Handling, Floating-Point Values on Special Numbers, MPFR Basics
 @comment  node-name,  next,  previous,  up
@@ -1581,7 +1585,8 @@
 @deftypefunx int mpfr_add_z (mpfr_t @var{rop}, mpfr_t @var{op1}, mpz_t @var{op2}, mpfr_rnd_t @var{rnd})
 @deftypefunx int mpfr_add_q (mpfr_t @var{rop}, mpfr_t @var{op1}, mpq_t @var{op2}, mpfr_rnd_t @var{rnd})
 Set @var{rop} to @math{@var{op1} + @var{op2}} rounded in the direction
-@var{rnd}. For types having no signed zero, it is considered unsigned
+@var{rnd}.  The IEEE-754 rules are used, in particular for signed zeros.
+But for types having no signed zeros, 0 is considered unsigned
 (i.e., (+0) + 0 = (+0) and (@minus{}0) + 0 = (@minus{}0)).
 The @code{mpfr_add_d} function assumes that the radix of the @code{double} type
 is a power of 2, with a precision at most that declared by the C implementation
@@ -1599,7 +1604,8 @@
 @deftypefunx int mpfr_sub_z (mpfr_t @var{rop}, mpfr_t @var{op1}, mpz_t @var{op2}, mpfr_rnd_t @var{rnd})
 @deftypefunx int mpfr_sub_q (mpfr_t @var{rop}, mpfr_t @var{op1}, mpq_t @var{op2}, mpfr_rnd_t @var{rnd})
 Set @var{rop} to @math{@var{op1} - @var{op2}} rounded in the direction
-@var{rnd}. For types having no signed zero, it is considered unsigned
+@var{rnd}.  The IEEE-754 rules are used, in particular for signed zeros.
+But for types having no signed zeros, 0 is considered unsigned
 (i.e., (+0) @minus{} 0 = (+0), (@minus{}0) @minus{} 0 = (@minus{}0),
 0 @minus{} (+0) = (@minus{}0) and 0 @minus{} (@minus{}0) = (+0)).
 The same restrictions than for @code{mpfr_add_d} apply to @code{mpfr_d_sub}
@@ -1615,7 +1621,7 @@
 Set @var{rop} to @math{@var{op1} @GMPtimes{} @var{op2}} rounded in the
 direction @var{rnd}.
 When a result is zero, its sign is the product of the signs of the operands
-(for types having no signed zero, it is considered positive).
+(for types having no signed zeros, 0 is considered positive).
 The same restrictions than for @code{mpfr_add_d} apply to @code{mpfr_mul_d}.
 @end deftypefun
 
@@ -1635,7 +1641,7 @@
 @deftypefunx int mpfr_div_q (mpfr_t @var{rop}, mpfr_t @var{op1}, mpq_t @var{op2}, mpfr_rnd_t @var{rnd})
 Set @var{rop} to @math{@var{op1}/@var{op2}} rounded in the direction @var{rnd}.
 When a result is zero, its sign is the product of the signs of the operands
-(for types having no signed zero, it is considered positive).
+(for types having no signed zeros, 0 is considered positive).
 The same restrictions than for @code{mpfr_add_d} apply to @code{mpfr_d_div}
 and @code{mpfr_div_d}.
 @end deftypefun
@@ -1643,15 +1649,18 @@
 @deftypefun int mpfr_sqrt (mpfr_t @var{rop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
 @deftypefunx int mpfr_sqrt_ui (mpfr_t @var{rop}, unsigned long int @var{op}, mpfr_rnd_t @var{rnd})
 Set @var{rop} to @m{\sqrt{@var{op}}, the square root of @var{op}}
-rounded in the direction @var{rnd} (set @var{rop} to @minus{}0 if @var{op} is
-@minus{}0, to be consistent with the IEEE 754 standard).
+rounded in the direction @var{rnd}.  Set @var{rop} to @minus{}0 if
+@var{op} is @minus{}0, to be consistent with the IEEE 754 standard.
 Set @var{rop} to NaN if @var{op} is negative.
 @end deftypefun
 
 @deftypefun int mpfr_rec_sqrt (mpfr_t @var{rop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
 Set @var{rop} to @m{1/\sqrt{@var{op}}, the reciprocal square root of @var{op}}
-rounded in the direction @var{rnd}. Set @var{rop} to +Inf if @var{op} is
-@pom{}0, +0 if @var{op} is +Inf, and NaN if @var{op} is negative.
+rounded in the direction @var{rnd}.  Set @var{rop} to +Inf if @var{op} is
+@pom{}0, +0 if @var{op} is +Inf, and NaN if @var{op} is negative.  Warning!
+Therefore the result on @minus{}0 is different from the one of the rSqrt
+function recommended by the IEEE 754-2008 standard (Section 9.2.1), which
+is @minus{}Inf instead of +Inf.
 @end deftypefun
 
 @deftypefun int mpfr_cbrt (mpfr_t @var{rop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
@@ -1832,7 +1841,9 @@
 @m{\log_2 @var{op}, log2(@var{op})} or
 @m{\log_{10} @var{op}, log10(@var{op})}, respectively,
 rounded in the direction @var{rnd}.
-Set @var{rop} to @minus{}Inf if @var{op} is @minus{}0
+Set @var{rop} to +0 if @var{op} is 1 (in all rounding modes),
+for consistency with the ISO C99 and IEEE 754-2008 standards.
+Set @var{rop} to @minus{}Inf if @var{op} is @pom{}0
 (i.e., the sign of the zero has no influence on the result).
 @end deftypefun
 
@@ -2003,8 +2014,11 @@
 @deftypefun int mpfr_lngamma (mpfr_t @var{rop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
 Set @var{rop} to the value of the logarithm of the Gamma function on @var{op},
 rounded in the direction @var{rnd}.
-When @math{@minus{}2@var{k}@minus{}1 @le{} @var{op} @le{} @minus{}2@var{k}},
-@var{k} being a non-negative integer, @var{rop} is set to NaN.
+When @var{op} is 1 or 2, set @var{rop} to +0 (in all rounding modes).
+When @var{op} is an infinity or a nonpositive integer, set @var{rop} to +Inf,
+following the general rules on special values.
+When @math{@minus{}2@var{k}@minus{}1 < @var{op} < @minus{}2@var{k}},
+@var{k} being a nonnegative integer, set @var{rop} to NaN@.
 See also @code{mpfr_lgamma}.
 @end deftypefun
 
@@ -2012,10 +2026,11 @@
 Set @var{rop} to the value of the logarithm of the absolute value of the
 Gamma function on @var{op}, rounded in the direction @var{rnd}. The sign
 (1 or @minus{}1) of Gamma(@var{op}) is returned in the object pointed to
-by @var{signp}. When @var{op} is an infinity or a non-positive integer, set
-@var{rop} to +Inf. When @var{op} is NaN, @minus{}Inf or a negative integer,
-*@var{signp} is undefined, and when @var{op} is @pom{}0, *@var{signp} is
-the sign of the zero.
+by @var{signp}.
+When @var{op} is 1 or 2, set @var{rop} to +0 (in all rounding modes).
+When @var{op} is an infinity or a nonpositive integer, set @var{rop} to +Inf.
+When @var{op} is NaN, @minus{}Inf or a negative integer, *@var{signp} is
+undefined, and when @var{op} is @pom{}0, *@var{signp} is the sign of the zero.
 @end deftypefun
 
 @deftypefun int mpfr_digamma (mpfr_t @var{rop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
@@ -2064,7 +2079,10 @@
 @deftypefunx int mpfr_fms (mpfr_t @var{rop}, mpfr_t @var{op1}, mpfr_t @var{op2}, mpfr_t @var{op3}, mpfr_rnd_t @var{rnd})
 Set @var{rop} to @math{(@var{op1} @GMPtimes{} @var{op2}) + @var{op3}}
 (resp.@: @math{(@var{op1} @GMPtimes{} @var{op2}) - @var{op3}})
-rounded in the direction @var{rnd}.
+rounded in the direction @var{rnd}.  Concerning special values (signed zeros,
+infinities, NaN), these functions behave like a multiplication followed by a
+separate addition or subtraction.  That is, the fused operation matters only
+for rounding.
 @end deftypefun
 
 @deftypefun int mpfr_agm (mpfr_t @var{rop}, mpfr_t @var{op1}, mpfr_t @var{op2}, mpfr_rnd_t @var{rnd})
@@ -2089,8 +2107,8 @@
 i.e., $\sqrt{x^2+y^2}$,
 @end tex
 rounded in the direction @var{rnd}.
-Special values are handled as described in Section F.9.4.3 of
-the ISO C99 and IEEE 754-2008 standards:
+Special values are handled as described in the ISO C99 (Section F.9.4.3)
+and IEEE 754-2008 (Section 9.2.1) standards:
 If @var{x} or @var{y} is an infinity, then +Inf is returned in @var{rop},
 even if the other number is NaN.
 @end deftypefun
--- a/PATCHES	2015-06-19 21:55:09.000000000 +0200
+++ b/PATCHES	2016-02-16 14:55:39.870392274 +0100
@@ -0,0 +1,12 @@
+root
+fits
+can_round
+si-ops
+sqrt
+zeta
+jn
+divhigh-basecase
+frexp
+muldiv-2exp-underflow
+muldiv-2exp-overflow
+lngamma-and-doc
--- a/src/div_2si.c	2015-06-19 21:55:10.000000000 +0200
+++ b/src/div_2si.c	2016-02-16 14:55:39.816391223 +0100
@@ -45,11 +45,12 @@
           if (rnd_mode == MPFR_RNDN &&
               (__gmpfr_emin > MPFR_EMAX_MAX - (n - 1) ||
                exp < __gmpfr_emin + (n - 1) ||
-               (inexact >= 0 && mpfr_powerof2_raw (y))))
+               ((MPFR_IS_NEG (y) ? inexact <= 0 : inexact >= 0) &&
+                mpfr_powerof2_raw (y))))
             rnd_mode = MPFR_RNDZ;
           return mpfr_underflow (y, rnd_mode, MPFR_SIGN(y));
         }
-      else if (MPFR_UNLIKELY(n < 0 && (__gmpfr_emax < MPFR_EMIN_MIN - n ||
+      else if (MPFR_UNLIKELY(n <= 0 && (__gmpfr_emax < MPFR_EMIN_MIN - n ||
                                        exp > __gmpfr_emax + n)) )
         return mpfr_overflow (y, rnd_mode, MPFR_SIGN(y));
 
--- a/src/div_2ui.c	2015-06-19 21:55:10.000000000 +0200
+++ b/src/div_2ui.c	2016-02-16 14:55:39.816391223 +0100
@@ -32,7 +32,7 @@
      rnd_mode),
     ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec(y), mpfr_log_prec, y, inexact));
 
-  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
+  if (MPFR_UNLIKELY (n == 0 || MPFR_IS_SINGULAR (x)))
     return mpfr_set (y, x, rnd_mode);
   else
     {
@@ -44,7 +44,9 @@
       if (MPFR_UNLIKELY (n >= diffexp))  /* exp - n <= emin - 1 */
         {
           if (rnd_mode == MPFR_RNDN &&
-              (n > diffexp || (inexact >= 0 && mpfr_powerof2_raw (y))))
+              (n > diffexp ||
+               ((MPFR_IS_NEG (y) ? inexact <= 0 : inexact >= 0) &&
+                mpfr_powerof2_raw (y))))
             rnd_mode = MPFR_RNDZ;
           return mpfr_underflow (y, rnd_mode, MPFR_SIGN (y));
         }
--- a/src/div.c	2015-06-19 21:55:09.000000000 +0200
+++ b/src/div.c	2016-02-16 14:55:39.857392021 +0100
@@ -310,24 +310,23 @@
 
       qp = MPFR_TMP_LIMBS_ALLOC (n);
       qh = mpfr_divhigh_n (qp, ap, bp, n);
+      MPFR_ASSERTD (qh == 0 || qh == 1);
       /* in all cases, the error is at most (2n+2) ulps on qh*B^n+{qp,n},
          cf algorithms.tex */
 
       p = n * GMP_NUMB_BITS - MPFR_INT_CEIL_LOG2 (2 * n + 2);
-      /* if qh is 1, then we need only PREC(q)-1 bits of {qp,n},
-         if rnd=RNDN, we need to be able to round with a directed rounding
-            and one more bit */
+      /* If rnd=RNDN, we need to be able to round with a directed rounding
+         and one more bit. */
+      if (qh == 1)
+        {
+          mpn_rshift (qp, qp, n, 1);
+          qp[n - 1] |= MPFR_LIMB_HIGHBIT;
+        }
       if (MPFR_LIKELY (mpfr_round_p (qp, n, p,
-                                 MPFR_PREC(q) + (rnd_mode == MPFR_RNDN) - qh)))
+                                     MPFR_PREC(q) + (rnd_mode == MPFR_RNDN))))
         {
           /* we can round correctly whatever the rounding mode */
-          if (qh == 0)
-            MPN_COPY (q0p, qp + 1, q0size);
-          else
-            {
-              mpn_rshift (q0p, qp + 1, q0size, 1);
-              q0p[q0size - 1] ^= MPFR_LIMB_HIGHBIT;
-            }
+          MPN_COPY (q0p, qp + 1, q0size);
           q0p[0] &= ~MPFR_LIMB_MASK(sh); /* put to zero low sh bits */
 
           if (rnd_mode == MPFR_RNDN) /* round to nearest */
@@ -335,15 +334,10 @@
               /* we know we can round, thus we are never in the even rule case:
                  if the round bit is 0, we truncate
                  if the round bit is 1, we add 1 */
-              if (qh == 0)
-                {
-                  if (sh > 0)
-                    round_bit = (qp[1] >> (sh - 1)) & 1;
-                  else
-                    round_bit = qp[0] >> (GMP_NUMB_BITS - 1);
-                }
-              else /* qh = 1 */
-                round_bit = (qp[1] >> sh) & 1;
+              if (sh > 0)
+                round_bit = (qp[1] >> (sh - 1)) & 1;
+              else
+                round_bit = qp[0] >> (GMP_NUMB_BITS - 1);
               if (round_bit == 0)
                 {
                   inex = -1;
--- a/src/div_ui.c	2015-06-19 21:55:10.000000000 +0200
+++ b/src/div_ui.c	2016-02-16 14:55:39.851391904 +0100
@@ -274,7 +274,8 @@
     res = mpfr_div_ui (y, x, u, rnd_mode);
   else
     {
-      res = -mpfr_div_ui (y, x, -u, MPFR_INVERT_RND (rnd_mode));
+      res = - mpfr_div_ui (y, x, - (unsigned long) u,
+                           MPFR_INVERT_RND (rnd_mode));
       MPFR_CHANGE_SIGN (y);
     }
   return res;
--- a/src/fits_intmax.c	2015-06-19 21:55:09.000000000 +0200
+++ b/src/fits_intmax.c	2016-02-16 14:55:39.864392157 +0100
@@ -33,6 +33,7 @@
 int
 mpfr_fits_intmax_p (mpfr_srcptr f, mpfr_rnd_t rnd)
 {
+  unsigned int saved_flags;
   mpfr_exp_t e;
   int prec;
   mpfr_t x, y;
@@ -85,6 +86,7 @@
   MPFR_ASSERTD (e == prec);
 
   /* hard case: first round to prec bits, then check */
+  saved_flags = __gmpfr_flags;
   mpfr_init2 (x, prec);
   mpfr_set (x, f, rnd);
 
@@ -97,10 +99,16 @@
     }
   else
     {
-      res = MPFR_GET_EXP (x) == e;
+      /* Warning! Due to the rounding, x can be an infinity. Here we use
+         the fact that singular numbers have a special exponent field,
+         thus well-defined and different from e, in which case this means
+         that the number does not fit. That's why we use MPFR_EXP, not
+         MPFR_GET_EXP. */
+      res = MPFR_EXP (x) == e;
     }
 
   mpfr_clear (x);
+  __gmpfr_flags = saved_flags;
   return res;
 }
 
--- a/src/fits_s.h	2015-06-19 21:55:10.000000000 +0200
+++ b/src/fits_s.h	2016-02-16 14:55:39.864392157 +0100
@@ -29,6 +29,7 @@
 int
 FUNCTION (mpfr_srcptr f, mpfr_rnd_t rnd)
 {
+  unsigned int saved_flags;
   mpfr_exp_t e;
   int prec;
   mpfr_t x;
@@ -81,9 +82,16 @@
   MPFR_ASSERTD (e == prec);
 
   /* hard case: first round to prec bits, then check */
+  saved_flags = __gmpfr_flags;
   mpfr_init2 (x, prec);
   mpfr_set (x, f, rnd);
-  res = neg ? (mpfr_cmp_si (x, MINIMUM) >= 0) : (MPFR_GET_EXP (x) == e);
+  /* Warning! Due to the rounding, x can be an infinity. Here we use
+     the fact that singular numbers have a special exponent field,
+     thus well-defined and different from e, in which case this means
+     that the number does not fit. That's why we use MPFR_EXP, not
+     MPFR_GET_EXP. */
+  res = neg ? (mpfr_cmp_si (x, MINIMUM) >= 0) : (MPFR_EXP (x) == e);
   mpfr_clear (x);
+  __gmpfr_flags = saved_flags;
   return res;
 }
--- a/src/fits_u.h	2015-06-19 21:55:09.000000000 +0200
+++ b/src/fits_u.h	2016-02-16 14:55:39.866392196 +0100
@@ -25,6 +25,7 @@
 int
 FUNCTION (mpfr_srcptr f, mpfr_rnd_t rnd)
 {
+  unsigned int saved_flags;
   mpfr_exp_t e;
   int prec;
   TYPE s;
@@ -62,9 +63,16 @@
   MPFR_ASSERTD (e == prec);
 
   /* hard case: first round to prec bits, then check */
+  saved_flags = __gmpfr_flags;
   mpfr_init2 (x, prec);
   mpfr_set (x, f, rnd);
-  res = MPFR_GET_EXP (x) == e;
+  /* Warning! Due to the rounding, x can be an infinity. Here we use
+     the fact that singular numbers have a special exponent field,
+     thus well-defined and different from e, in which case this means
+     that the number does not fit. That's why we use MPFR_EXP, not
+     MPFR_GET_EXP. */
+  res = MPFR_EXP (x) == e;
   mpfr_clear (x);
+  __gmpfr_flags = saved_flags;
   return res;
 }
--- a/src/frexp.c	2015-06-19 21:55:09.000000000 +0200
+++ b/src/frexp.c	2016-02-16 14:55:39.820391301 +0100
@@ -26,6 +26,13 @@
 mpfr_frexp (mpfr_exp_t *exp, mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
 {
   int inex;
+  unsigned int saved_flags = __gmpfr_flags;
+  MPFR_BLOCK_DECL (flags);
+
+  MPFR_LOG_FUNC
+    (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd),
+     ("y[%Pu]=%.*Rg exp=%" MPFR_EXP_FSPEC "d inex=%d", mpfr_get_prec (y),
+      mpfr_log_prec, y, (mpfr_eexp_t) *exp, inex));
 
   if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x)))
     {
@@ -49,8 +56,32 @@
         }
     }
 
-  inex = mpfr_set (y, x, rnd);
+  MPFR_BLOCK (flags, inex = mpfr_set (y, x, rnd));
+  __gmpfr_flags = saved_flags;
+
+  /* Possible overflow due to the rounding, no possible underflow. */
+
+  if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
+    {
+      int inex2;
+
+      /* An overflow here means that the exponent of y would be larger than
+         the one of x, thus x would be rounded to the next power of 2, and
+         the returned y should be 1/2 in absolute value, rounded (i.e. with
+         possible underflow or overflow). This also implies that x and y are
+         different objects, so that the exponent of x has not been lost. */
+      MPFR_LOG_MSG (("Internal overflow\n", 0));
+      MPFR_ASSERTD (x != y);
+      *exp = MPFR_GET_EXP (x) + 1;
+      inex2 = mpfr_set_si_2exp (y, MPFR_INT_SIGN (x), -1, rnd);
+      MPFR_LOG_MSG (("inex=%d inex2=%d\n", inex, inex2));
+      if (inex2 != 0)
+        inex = inex2;
+      MPFR_RET (inex);
+    }
+
   *exp = MPFR_GET_EXP (y);
-  MPFR_SET_EXP (y, 0);
+  /* Do not use MPFR_SET_EXP because the range has not been checked yet. */
+  MPFR_EXP (y) = 0;
   return mpfr_check_range (y, inex, rnd);
 }
--- a/src/jyn_asympt.c	2015-06-19 21:55:09.000000000 +0200
+++ b/src/jyn_asympt.c	2016-02-16 14:55:39.829391476 +0100
@@ -253,9 +253,9 @@
         break;
       if (diverge != 0)
         {
-          mpfr_set (c, z, r); /* will force inex=0 below, which means the
-                               asymptotic expansion failed */
-          break;
+          MPFR_ZIV_FREE (loop);
+          mpfr_clear (c);
+          return 0; /* means that the asymptotic expansion failed */
         }
       MPFR_ZIV_NEXT (loop, w);
     }
--- a/src/lngamma.c	2015-06-19 21:55:10.000000000 +0200
+++ b/src/lngamma.c	2016-02-16 14:55:39.807391048 +0100
@@ -603,16 +603,17 @@
       mpfr_get_prec (y), mpfr_log_prec, y, inex));
 
   /* special cases */
-  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
+  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x) ||
+                     (MPFR_IS_NEG (x) && mpfr_integer_p (x))))
     {
-      if (MPFR_IS_NAN (x) || MPFR_IS_NEG (x))
+      if (MPFR_IS_NAN (x))
         {
           MPFR_SET_NAN (y);
           MPFR_RET_NAN;
         }
-      else /* lngamma(+Inf) = lngamma(+0) = +Inf */
+      else /* lngamma(+/-Inf) = lngamma(nonpositive integer) = +Inf */
         {
-          if (MPFR_IS_ZERO (x))
+          if (!MPFR_IS_INF (x))
             mpfr_set_divby0 ();
           MPFR_SET_INF (y);
           MPFR_SET_POS (y);
@@ -620,8 +621,8 @@
         }
     }
 
-  /* if x < 0 and -2k-1 <= x <= -2k, then lngamma(x) = NaN */
-  if (MPFR_IS_NEG (x) && (unit_bit (x) == 0 || mpfr_integer_p (x)))
+  /* if -2k-1 < x < -2k <= 0, then lngamma(x) = NaN */
+  if (MPFR_IS_NEG (x) && unit_bit (x) == 0)
     {
       MPFR_SET_NAN (y);
       MPFR_RET_NAN;
--- a/src/mpfr.h	2015-06-19 21:55:10.000000000 +0200
+++ b/src/mpfr.h	2016-02-16 14:55:39.872392313 +0100
@@ -27,7 +27,7 @@
 #define MPFR_VERSION_MAJOR 3
 #define MPFR_VERSION_MINOR 1
 #define MPFR_VERSION_PATCHLEVEL 3
-#define MPFR_VERSION_STRING "3.1.3"
+#define MPFR_VERSION_STRING "3.1.3-p12"
 
 /* Macros dealing with MPFR VERSION */
 #define MPFR_VERSION_NUM(a,b,c) (((a) << 16L) | ((b) << 8) | (c))
--- a/src/mul_2si.c	2015-06-19 21:55:10.000000000 +0200
+++ b/src/mul_2si.c	2016-02-16 14:55:39.817391242 +0100
@@ -39,7 +39,7 @@
     {
       mpfr_exp_t exp = MPFR_GET_EXP (x);
       MPFR_SETRAW (inexact, y, x, exp, rnd_mode);
-      if (MPFR_UNLIKELY( n > 0 && (__gmpfr_emax < MPFR_EMIN_MIN + n ||
+      if (MPFR_UNLIKELY(n >= 0 && (__gmpfr_emax < MPFR_EMIN_MIN + n ||
                                    exp > __gmpfr_emax - n)))
         return mpfr_overflow (y, rnd_mode, MPFR_SIGN(y));
       else if (MPFR_UNLIKELY(n < 0 && (__gmpfr_emin > MPFR_EMAX_MAX + n ||
@@ -48,7 +48,8 @@
           if (rnd_mode == MPFR_RNDN &&
               (__gmpfr_emin > MPFR_EMAX_MAX + (n + 1) ||
                exp < __gmpfr_emin - (n + 1) ||
-               (inexact >= 0 && mpfr_powerof2_raw (y))))
+               ((MPFR_IS_NEG (y) ? inexact <= 0 : inexact >= 0) &&
+                mpfr_powerof2_raw (y))))
             rnd_mode = MPFR_RNDZ;
           return mpfr_underflow (y, rnd_mode, MPFR_SIGN(y));
         }
--- a/src/mulders.c	2015-06-19 21:55:10.000000000 +0200
+++ b/src/mulders.c	2016-02-16 14:55:39.827391437 +0100
@@ -236,9 +236,10 @@
          that in addition to the limb np[n-1] to reduce, we have at least 2
          extra limbs, thus accessing np[n-3] is valid. */
 
-      /* warning: we can have np[n-1]=d1 and np[n-2]=d0, but since {np,n} < D,
-         the largest possible partial quotient is B-1 */
-      if (MPFR_UNLIKELY(np[n - 1] == d1 && np[n - 2] == d0))
+      /* Warning: we can have np[n-1]>d1 or (np[n-1]=d1 and np[n-2]>=d0) here,
+         since we truncate the divisor at each step, but since {np,n} < D
+         originally, the largest possible partial quotient is B-1. */
+      if (MPFR_UNLIKELY(np[n-1] > d1 || (np[n-1] == d1 && np[n-2] >= d0)))
         q2 = ~ (mp_limb_t) 0;
       else
         udiv_qr_3by2 (q2, q1, q0, np[n - 1], np[n - 2], np[n - 3],
--- a/src/mul_ui.c	2015-06-19 21:55:10.000000000 +0200
+++ b/src/mul_ui.c	2016-02-16 14:55:39.854391963 +0100
@@ -126,7 +126,8 @@
     res = mpfr_mul_ui (y, x, u, rnd_mode);
   else
     {
-      res = -mpfr_mul_ui (y, x, -u, MPFR_INVERT_RND (rnd_mode));
+      res = - mpfr_mul_ui (y, x, - (unsigned long) u,
+                           MPFR_INVERT_RND (rnd_mode));
       MPFR_CHANGE_SIGN (y);
     }
   return res;
--- a/src/root.c	2015-06-19 21:55:10.000000000 +0200
+++ b/src/root.c	2016-02-16 14:55:39.873392333 +0100
@@ -23,13 +23,15 @@
 #define MPFR_NEED_LONGLONG_H
 #include "mpfr-impl.h"
 
- /* The computation of y = x^(1/k) is done as follows:
+ /* The computation of y = x^(1/k) is done as follows, except for large
+    values of k, for which this would be inefficient or yield internal
+    integer overflows:
 
     Let x = sign * m * 2^(k*e) where m is an integer
 
     with 2^(k*(n-1)) <= m < 2^(k*n) where n = PREC(y)
 
-    and m = s^k + r where 0 <= r and m < (s+1)^k
+    and m = s^k + t where 0 <= t and m < (s+1)^k
 
     we want that s has n bits i.e. s >= 2^(n-1), or m >= 2^(k*(n-1))
     i.e. m must have at least k*(n-1)+1 bits
@@ -38,11 +40,15 @@
     x^(1/k) = s * 2^e or (s+1) * 2^e according to the rounding mode.
  */
 
+static int
+mpfr_root_aux (mpfr_ptr y, mpfr_srcptr x, unsigned long k,
+               mpfr_rnd_t rnd_mode);
+
 int
 mpfr_root (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode)
 {
   mpz_t m;
-  mpfr_exp_t e, r, sh;
+  mpfr_exp_t e, r, sh, f;
   mpfr_prec_t n, size_m, tmp;
   int inexact, negative;
   MPFR_SAVE_EXPO_DECL (expo);
@@ -55,50 +61,27 @@
 
   if (MPFR_UNLIKELY (k <= 1))
     {
-      if (k < 1) /* k==0 => y=x^(1/0)=x^(+Inf) */
-#if 0
-        /* For 0 <= x < 1 => +0.
-           For x = 1      => 1.
-           For x > 1,     => +Inf.
-           For x < 0      => NaN.
-        */
+      if (k == 0)
         {
-          if (MPFR_IS_NEG (x) && !MPFR_IS_ZERO (x))
-            {
-              MPFR_SET_NAN (y);
-              MPFR_RET_NAN;
-            }
-          inexact = mpfr_cmp (x, __gmpfr_one);
-          if (inexact == 0)
-            return mpfr_set_ui (y, 1, rnd_mode); /* 1 may be Out of Range */
-          else if (inexact < 0)
-            return mpfr_set_ui (y, 0, rnd_mode); /* 0+ */
-          else
-            {
-              mpfr_set_inf (y, 1);
-              return 0;
-            }
+          MPFR_SET_NAN (y);
+          MPFR_RET_NAN;
         }
-#endif
-      {
-        MPFR_SET_NAN (y);
-        MPFR_RET_NAN;
-      }
-      else /* y =x^(1/1)=x */
+      else /* y = x^(1/1) = x */
         return mpfr_set (y, x, rnd_mode);
     }
 
   /* Singular values */
-  else if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
+  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
     {
       if (MPFR_IS_NAN (x))
         {
           MPFR_SET_NAN (y); /* NaN^(1/k) = NaN */
           MPFR_RET_NAN;
         }
-      else if (MPFR_IS_INF (x)) /* +Inf^(1/k) = +Inf
-                                   -Inf^(1/k) = -Inf if k odd
-                                   -Inf^(1/k) = NaN if k even */
+
+      if (MPFR_IS_INF (x)) /* +Inf^(1/k) = +Inf
+                              -Inf^(1/k) = -Inf if k odd
+                              -Inf^(1/k) = NaN if k even */
         {
           if (MPFR_IS_NEG(x) && (k % 2 == 0))
             {
@@ -106,27 +89,31 @@
               MPFR_RET_NAN;
             }
           MPFR_SET_INF (y);
-          MPFR_SET_SAME_SIGN (y, x);
-          MPFR_RET (0);
         }
       else /* x is necessarily 0: (+0)^(1/k) = +0
                                   (-0)^(1/k) = -0 */
         {
           MPFR_ASSERTD (MPFR_IS_ZERO (x));
           MPFR_SET_ZERO (y);
-          MPFR_SET_SAME_SIGN (y, x);
-          MPFR_RET (0);
         }
+      MPFR_SET_SAME_SIGN (y, x);
+      MPFR_RET (0);
     }
 
   /* Returns NAN for x < 0 and k even */
-  else if (MPFR_IS_NEG (x) && (k % 2 == 0))
+  if (MPFR_UNLIKELY (MPFR_IS_NEG (x) && (k % 2 == 0)))
     {
       MPFR_SET_NAN (y);
       MPFR_RET_NAN;
     }
 
   /* General case */
+
+  /* For large k, use exp(log(x)/k). The threshold of 100 seems to be quite
+     good when the precision goes to infinity. */
+  if (k > 100)
+    return mpfr_root_aux (y, x, k, rnd_mode);
+
   MPFR_SAVE_EXPO_MARK (expo);
   mpz_init (m);
 
@@ -135,31 +122,24 @@
     mpz_neg (m, m);
   r = e % (mpfr_exp_t) k;
   if (r < 0)
-    r += k; /* now r = e (mod k) with 0 <= e < r */
+    r += k; /* now r = e (mod k) with 0 <= r < k */
+  MPFR_ASSERTD (0 <= r && r < k);
   /* x = (m*2^r) * 2^(e-r) where e-r is a multiple of k */
 
   MPFR_MPZ_SIZEINBASE2 (size_m, m);
   /* for rounding to nearest, we want the round bit to be in the root */
   n = MPFR_PREC (y) + (rnd_mode == MPFR_RNDN);
 
-  /* we now multiply m by 2^(r+k*sh) so that root(m,k) will give
-     exactly n bits: we want k*(n-1)+1 <= size_m + k*sh + r <= k*n
-     i.e. sh = floor ((kn-size_m-r)/k) */
-  if ((mpfr_exp_t) size_m + r > k * (mpfr_exp_t) n)
-    sh = 0; /* we already have too many bits */
+  /* we now multiply m by 2^sh so that root(m,k) will give
+     exactly n bits: we want k*(n-1)+1 <= size_m + sh <= k*n
+     i.e. sh = k*f + r with f = max(floor((k*n-size_m-r)/k),0) */
+  if ((mpfr_exp_t) size_m + r >= k * (mpfr_exp_t) n)
+    f = 0; /* we already have too many bits */
   else
-    sh = (k * (mpfr_exp_t) n - (mpfr_exp_t) size_m - r) / k;
-  sh = k * sh + r;
-  if (sh >= 0)
-    {
-      mpz_mul_2exp (m, m, sh);
-      e = e - sh;
-    }
-  else if (r > 0)
-    {
-      mpz_mul_2exp (m, m, r);
-      e = e - r;
-    }
+    f = (k * (mpfr_exp_t) n - (mpfr_exp_t) size_m - r) / k;
+  sh = k * f + r;
+  mpz_mul_2exp (m, m, sh);
+  e = e - sh;
 
   /* invariant: x = m*2^e, with e divisible by k */
 
@@ -203,3 +183,97 @@
   MPFR_SAVE_EXPO_FREE (expo);
   return mpfr_check_range (y, inexact, rnd_mode);
 }
+
+/* Compute y <- x^(1/k) using exp(log(x)/k).
+   Assume all special cases have been eliminated before.
+   In the extended exponent range, overflows/underflows are not possible.
+   Assume x > 0, or x < 0 and k odd.
+*/
+static int
+mpfr_root_aux (mpfr_ptr y, mpfr_srcptr x, unsigned long k, mpfr_rnd_t rnd_mode)
+{
+  int inexact, exact_root = 0;
+  mpfr_prec_t w; /* working precision */
+  mpfr_t absx, t;
+  MPFR_GROUP_DECL(group);
+  MPFR_TMP_DECL(marker);
+  MPFR_ZIV_DECL(loop);
+  MPFR_SAVE_EXPO_DECL (expo);
+
+  MPFR_TMP_INIT_ABS (absx, x);
+
+  MPFR_TMP_MARK(marker);
+  w = MPFR_PREC(y) + 10;
+  /* Take some guard bits to prepare for the 'expt' lost bits below.
+     If |x| < 2^k, then log|x| < k, thus taking log2(k) bits should be fine. */
+  if (MPFR_GET_EXP(x) > 0)
+    w += MPFR_INT_CEIL_LOG2 (MPFR_GET_EXP(x));
+  MPFR_GROUP_INIT_1(group, w, t);
+  MPFR_SAVE_EXPO_MARK (expo);
+  MPFR_ZIV_INIT (loop, w);
+  for (;;)
+    {
+      mpfr_exp_t expt;
+      unsigned int err;
+
+      mpfr_log (t, absx, MPFR_RNDN);
+      /* t = log|x| * (1 + theta) with |theta| <= 2^(-w) */
+      mpfr_div_ui (t, t, k, MPFR_RNDN);
+      expt = MPFR_GET_EXP (t);
+      /* t = log|x|/k * (1 + theta) + eps with |theta| <= 2^(-w)
+         and |eps| <= 1/2 ulp(t), thus the total error is bounded
+         by 1.5 * 2^(expt - w) */
+      mpfr_exp (t, t, MPFR_RNDN);
+      /* t = |x|^(1/k) * exp(tau) * (1 + theta1) with
+         |tau| <= 1.5 * 2^(expt - w) and |theta1| <= 2^(-w).
+         For |tau| <= 0.5 we have |exp(tau)-1| < 4/3*tau, thus
+         for w >= expt + 2 we have:
+         t = |x|^(1/k) * (1 + 2^(expt+2)*theta2) * (1 + theta1) with
+         |theta1|, |theta2| <= 2^(-w).
+         If expt+2 > 0, as long as w >= 1, we have:
+         t = |x|^(1/k) * (1 + 2^(expt+3)*theta3) with |theta3| < 2^(-w).
+         For expt+2 = 0, we have:
+         t = |x|^(1/k) * (1 + 2^2*theta3) with |theta3| < 2^(-w).
+         Finally for expt+2 < 0 we have:
+         t = |x|^(1/k) * (1 + 2*theta3) with |theta3| < 2^(-w).
+      */
+      err = (expt + 2 > 0) ? expt + 3
+        : (expt + 2 == 0) ? 2 : 1;
+      /* now t = |x|^(1/k) * (1 + 2^(err-w)) thus the error is at most
+         2^(EXP(t) - w + err) */
+      if (MPFR_LIKELY (MPFR_CAN_ROUND(t, w - err, MPFR_PREC(y), rnd_mode)))
+        break;
+
+      /* If we fail to round correctly, check for an exact result or a
+         midpoint result with MPFR_RNDN (regarded as hard-to-round in
+         all precisions in order to determine the ternary value). */
+      {
+        mpfr_t z, zk;
+
+        mpfr_init2 (z, MPFR_PREC(y) + (rnd_mode == MPFR_RNDN));
+        mpfr_init2 (zk, MPFR_PREC(x));
+        mpfr_set (z, t, MPFR_RNDN);
+        inexact = mpfr_pow_ui (zk, z, k, MPFR_RNDN);
+        exact_root = !inexact && mpfr_equal_p (zk, absx);
+        if (exact_root) /* z is the exact root, thus round z directly */
+          inexact = mpfr_set4 (y, z, rnd_mode, MPFR_SIGN (x));
+        mpfr_clear (zk);
+        mpfr_clear (z);
+        if (exact_root)
+          break;
+      }
+
+      MPFR_ZIV_NEXT (loop, w);
+      MPFR_GROUP_REPREC_1(group, w, t);
+    }
+  MPFR_ZIV_FREE (loop);
+
+  if (!exact_root)
+    inexact = mpfr_set4 (y, t, rnd_mode, MPFR_SIGN (x));
+
+  MPFR_GROUP_CLEAR(group);
+  MPFR_TMP_FREE(marker);
+  MPFR_SAVE_EXPO_FREE (expo);
+
+  return mpfr_check_range (y, inexact, rnd_mode);
+}
--- a/src/round_p.c	2015-06-19 21:55:09.000000000 +0200
+++ b/src/round_p.c	2016-02-16 14:55:39.860392080 +0100
@@ -31,7 +31,11 @@
 {
   int i1, i2;
 
+  MPFR_ASSERTN(bp[bn - 1] & MPFR_LIMB_HIGHBIT);
+
   i1 = mpfr_round_p_2 (bp, bn, err0, prec);
+
+  /* compare with mpfr_can_round_raw */
   i2 = mpfr_can_round_raw (bp, bn, MPFR_SIGN_POS, err0,
                            MPFR_RNDN, MPFR_RNDZ, prec);
   if (i1 != i2)
@@ -42,6 +46,7 @@
       gmp_fprintf (stderr, "%NX\n", bp, bn);
       MPFR_ASSERTN (0);
     }
+
   return i1;
 }
 # define mpfr_round_p mpfr_round_p_2
@@ -62,6 +67,8 @@
   mp_limb_t tmp, mask;
   int s;
 
+  MPFR_ASSERTD(bp[bn - 1] & MPFR_LIMB_HIGHBIT);
+
   err = (mpfr_prec_t) bn * GMP_NUMB_BITS;
   if (MPFR_UNLIKELY (err0 <= 0 || (mpfr_uexp_t) err0 <= prec || prec >= err))
     return 0;  /* can't round */
--- a/src/round_prec.c	2015-06-19 21:55:10.000000000 +0200
+++ b/src/round_prec.c	2016-02-16 14:55:39.860392080 +0100
@@ -141,24 +141,40 @@
 mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err0,
                     mpfr_rnd_t rnd1, mpfr_rnd_t rnd2, mpfr_prec_t prec)
 {
-  mpfr_prec_t err;
+  mpfr_prec_t err, prec0 = prec;
   mp_size_t k, k1, tn;
   int s, s1;
   mp_limb_t cc, cc2;
   mp_limb_t *tmp;
   MPFR_TMP_DECL(marker);
 
+  MPFR_ASSERTD(bp[bn - 1] & MPFR_LIMB_HIGHBIT);
+
   if (MPFR_UNLIKELY(err0 < 0 || (mpfr_uexp_t) err0 <= prec))
     return 0;  /* can't round */
-  else if (MPFR_UNLIKELY (prec > (mpfr_prec_t) bn * GMP_NUMB_BITS))
-    { /* then ulp(b) < precision < error */
-      return rnd2 == MPFR_RNDN && (mpfr_uexp_t) err0 - 2 >= prec;
-      /* can round only in rounding to the nearest and err0 >= prec + 2 */
-    }
 
   MPFR_ASSERT_SIGN(neg);
   neg = MPFR_IS_NEG_SIGN(neg);
 
+  /* Transform RNDD and RNDU to Zero / Away */
+  MPFR_ASSERTD((neg == 0) || (neg == 1));
+  if (rnd1 != MPFR_RNDN)
+    rnd1 = MPFR_IS_LIKE_RNDZ(rnd1, neg) ? MPFR_RNDZ : MPFR_RNDA;
+  if (rnd2 != MPFR_RNDN)
+    rnd2 = MPFR_IS_LIKE_RNDZ(rnd2, neg) ? MPFR_RNDZ : MPFR_RNDA;
+
+  if (MPFR_UNLIKELY (prec > (mpfr_prec_t) bn * GMP_NUMB_BITS))
+    { /* Then prec < PREC(b): we can round:
+         (i) in rounding to the nearest iff err0 >= prec + 2
+         (ii) in directed rounding mode iff rnd1 is compatible with rnd2
+              and err0 >= prec + 1, unless b = 2^k and rnd1=rnd2=RNDA in
+              which case we need err0 >= prec + 2. */
+      if (rnd2 == MPFR_RNDN)
+        return (mpfr_uexp_t) err0 - 2 >= prec;
+      else
+        return (rnd1 == rnd2) && (mpfr_uexp_t) err0 - 2 >= prec;
+    }
+
   /* if the error is smaller than ulp(b), then anyway it will propagate
      up to ulp(b) */
   err = ((mpfr_uexp_t) err0 > (mpfr_prec_t) bn * GMP_NUMB_BITS) ?
@@ -168,19 +184,25 @@
   k = (err - 1) / GMP_NUMB_BITS;
   MPFR_UNSIGNED_MINUS_MODULO(s, err);
   /* the error corresponds to bit s in limb k, the most significant limb
-     being limb 0 */
+     being limb 0; in memory, limb k is bp[bn-1-k]. */
 
   k1 = (prec - 1) / GMP_NUMB_BITS;
   MPFR_UNSIGNED_MINUS_MODULO(s1, prec);
-  /* the last significant bit is bit s1 in limb k1 */
+  /* the least significant bit is bit s1 in limb k1 */
 
-  /* don't need to consider the k1 most significant limbs */
+  /* We don't need to consider the k1 most significant limbs.
+     They will be considered later only to detect when subtracting
+     the error bound yields a change of binade.
+     Warning! The number with updated bn may no longer be normalized. */
   k -= k1;
   bn -= k1;
   prec -= (mpfr_prec_t) k1 * GMP_NUMB_BITS;
 
-  /* if when adding or subtracting (1 << s) in bp[bn-1-k], it does not
-     change bp[bn-1] >> s1, then we can round */
+  /* We can decide of the correct rounding if rnd2(b-eps) and rnd2(b+eps)
+     give the same result to the target precision 'prec', i.e., if when
+     adding or subtracting (1 << s) in bp[bn-1-k], it does not change the
+     rounding in direction 'rnd2' at ulp-position bp[bn-1] >> s1, taking also
+     into account the possible change of binade. */
   MPFR_TMP_MARK(marker);
   tn = bn;
   k++; /* since we work with k+1 everywhere */
@@ -190,11 +212,6 @@
 
   MPFR_ASSERTD (k > 0);
 
-  /* Transform RNDD and RNDU to Zero / Away */
-  MPFR_ASSERTD((neg == 0) || (neg ==1));
-  if (MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd1, neg))
-    rnd1 = MPFR_RNDZ;
-
   switch (rnd1)
     {
     case MPFR_RNDZ:
@@ -203,33 +220,54 @@
       /* mpfr_round_raw2 returns 1 if one should add 1 at ulp(b,prec),
          and 0 otherwise */
       cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec);
-      /* cc is the new value of bit s1 in bp[bn-1] */
+      /* cc is the new value of bit s1 in bp[bn-1] after rounding 'rnd2' */
+
       /* now round b + 2^(MPFR_EXP(b)-err) */
-      cc2 = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
+      mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
+      /* if there was a carry here, then necessarily bit s1 of bp[bn-1]
+         changed, thus we surely cannot round for directed rounding, but this
+         will be detected below, with cc2 != cc */
       break;
     case MPFR_RNDN:
       /* Round to nearest */
-       /* first round b+2^(MPFR_EXP(b)-err) */
-      cc = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
+
+      /* first round b+2^(MPFR_EXP(b)-err) */
+      mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
+      /* same remark as above in case a carry occurs in mpn_add_1() */
       cc = (tmp[bn - 1] >> s1) & 1; /* gives 0 when cc=1 */
       cc ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec);
+      /* cc is the new value of bit s1 in bp[bn-1]+eps after rounding 'rnd2' */
+
+    subtract_eps:
       /* now round b-2^(MPFR_EXP(b)-err) */
       cc2 = mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
+      /* propagate the potential borrow up to the most significant limb
+         (it cannot propagate further since the most significant limb is
+         at least MPFR_LIMB_HIGHBIT) */
+      for (tn = 0; tn + 1 < k1 && (cc2 != 0); tn ++)
+        cc2 = bp[bn + tn] == 0;
+      /* We have an exponent decrease when either:
+           (i) k1 = 0 and tmp[bn-1] < MPFR_LIMB_HIGHBIT
+           (ii) k1 > 0 and cc <> 0 and bp[bn + tn] = MPFR_LIMB_HIGHBIT
+                (then necessarily tn = k1-1).
+         Then for directed rounding we cannot round,
+         and for rounding to nearest we cannot round when err = prec + 1.
+      */
+      if (((k1 == 0 && tmp[bn - 1] < MPFR_LIMB_HIGHBIT) ||
+           (k1 != 0 && cc2 != 0 && bp[bn + tn] == MPFR_LIMB_HIGHBIT)) &&
+          (rnd2 != MPFR_RNDN || err0 == prec0 + 1))
+        {
+          MPFR_TMP_FREE(marker);
+          return 0;
+        }
       break;
     default:
       /* Round away */
       cc = (bp[bn - 1] >> s1) & 1;
       cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec);
-      /* now round b +/- 2^(MPFR_EXP(b)-err) */
-      cc2 = mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
-      break;
-    }
+      /* cc is the new value of bit s1 in bp[bn-1]+eps after rounding 'rnd2' */
 
-  /* if cc2 is 1, then a carry or borrow propagates to the next limb */
-  if (cc2 && cc)
-    {
-      MPFR_TMP_FREE(marker);
-      return 0;
+      goto subtract_eps;
     }
 
   cc2 = (tmp[bn - 1] >> s1) & 1;
--- a/src/si_op.c	2015-06-19 21:55:10.000000000 +0200
+++ b/src/si_op.c	2016-02-16 14:55:39.854391963 +0100
@@ -30,7 +30,7 @@
   if (u >= 0)
     return mpfr_add_ui (y, x, u, rnd_mode);
   else
-    return mpfr_sub_ui (y, x, -u, rnd_mode);
+    return mpfr_sub_ui (y, x, - (unsigned long) u, rnd_mode);
 }
 
 int
@@ -39,7 +39,7 @@
   if (u >= 0)
     return mpfr_sub_ui (y, x, u, rnd_mode);
   else
-    return mpfr_add_ui (y, x, -u, rnd_mode);
+    return mpfr_add_ui (y, x, - (unsigned long) u, rnd_mode);
 }
 
 int
@@ -49,9 +49,9 @@
     return mpfr_ui_sub (y, u, x, rnd_mode);
   else
     {
-    int res = -mpfr_add_ui (y, x, -u, MPFR_INVERT_RND (rnd_mode));
-    MPFR_CHANGE_SIGN (y);
-    return res;
+      int res = - mpfr_add_ui (y, x, - (unsigned long) u,
+                               MPFR_INVERT_RND (rnd_mode));
+      MPFR_CHANGE_SIGN (y);
+      return res;
     }
 }
-
--- a/src/sqrt.c	2015-06-19 21:55:09.000000000 +0200
+++ b/src/sqrt.c	2016-02-16 14:55:39.836391612 +0100
@@ -211,10 +211,11 @@
       rsize --;
       sh = 0;
     }
+  /* now rsize = MPFR_LIMB_SIZE(r) */
   if (mpn_add_1 (rp0, rp, rsize, MPFR_LIMB_ONE << sh))
     {
       expr ++;
-      rp[rsize - 1] = MPFR_LIMB_HIGHBIT;
+      rp0[rsize - 1] = MPFR_LIMB_HIGHBIT;
     }
   goto end;
 
--- a/src/ui_div.c	2015-06-19 21:55:10.000000000 +0200
+++ b/src/ui_div.c	2016-02-16 14:55:39.854391963 +0100
@@ -106,7 +106,8 @@
     res = mpfr_ui_div (y, u, x, rnd_mode);
   else
     {
-      res = -mpfr_ui_div (y, -u, x, MPFR_INVERT_RND(rnd_mode));
+      res = - mpfr_ui_div (y, - (unsigned long) u, x,
+                           MPFR_INVERT_RND(rnd_mode));
       MPFR_CHANGE_SIGN (y);
     }
   return res;
--- a/src/version.c	2015-06-19 21:55:10.000000000 +0200
+++ b/src/version.c	2016-02-16 14:55:39.875392372 +0100
@@ -25,5 +25,5 @@
 const char *
 mpfr_get_version (void)
 {
-  return "3.1.3";
+  return "3.1.3-p12";
 }
--- a/src/zeta.c	2015-06-19 21:55:10.000000000 +0200
+++ b/src/zeta.c	2016-02-16 14:55:39.833391554 +0100
@@ -377,8 +377,8 @@
         }
     }
 
-  /* Check for case s= 1 before changing the exponent range */
-  if (mpfr_cmp (s, __gmpfr_one) ==0)
+  /* Check for case s=1 before changing the exponent range */
+  if (mpfr_cmp (s, __gmpfr_one) == 0)
     {
       MPFR_SET_INF (z);
       MPFR_SET_POS (z);
@@ -420,7 +420,7 @@
       MPFR_ZIV_INIT (loop, prec1);
       for (;;)
         {
-          mpfr_sub (s1, __gmpfr_one, s, MPFR_RNDN);/* s1 = 1-s */
+          mpfr_sub (s1, __gmpfr_one, s, MPFR_RNDN); /* s1 = 1-s */
           mpfr_zeta_pos (z_pre, s1, MPFR_RNDN);   /* zeta(1-s)  */
           mpfr_gamma (y, s1, MPFR_RNDN);          /* gamma(1-s) */
           if (MPFR_IS_INF (y)) /* Zeta(s) < 0 for -4k-2 < s < -4k,
@@ -432,17 +432,32 @@
               break;
             }
           mpfr_mul (z_pre, z_pre, y, MPFR_RNDN);  /* gamma(1-s)*zeta(1-s) */
-          mpfr_const_pi (p, MPFR_RNDD);
-          mpfr_mul (y, s, p, MPFR_RNDN);
-          mpfr_div_2ui (y, y, 1, MPFR_RNDN);      /* s*Pi/2 */
-          mpfr_sin (y, y, MPFR_RNDN);             /* sin(Pi*s/2) */
-          mpfr_mul (z_pre, z_pre, y, MPFR_RNDN);
+
+          mpfr_const_pi (p, MPFR_RNDD); /* p is Pi */
+
+          /* multiply z_pre by 2^s*Pi^(s-1) where p=Pi, s1=1-s */
           mpfr_mul_2ui (y, p, 1, MPFR_RNDN);      /* 2*Pi */
           mpfr_neg (s1, s1, MPFR_RNDN);           /* s-1 */
           mpfr_pow (y, y, s1, MPFR_RNDN);         /* (2*Pi)^(s-1) */
           mpfr_mul (z_pre, z_pre, y, MPFR_RNDN);
           mpfr_mul_2ui (z_pre, z_pre, 1, MPFR_RNDN);
 
+          /* multiply z_pre by sin(Pi*s/2) */
+          mpfr_mul (y, s, p, MPFR_RNDN);
+          mpfr_div_2ui (p, y, 1, MPFR_RNDN);      /* p = s*Pi/2 */
+          mpfr_sin (y, p, MPFR_RNDN);             /* y = sin(Pi*s/2) */
+          if (MPFR_GET_EXP(y) < 0) /* take account of cancellation in sin(p) */
+            {
+              mpfr_t t;
+              mpfr_init2 (t, prec1 - MPFR_GET_EXP(y));
+              mpfr_const_pi (t, MPFR_RNDD);
+              mpfr_mul (t, s, t, MPFR_RNDN);
+              mpfr_div_2ui (t, t, 1, MPFR_RNDN);
+              mpfr_sin (y, t, MPFR_RNDN);
+              mpfr_clear (t);
+            }
+          mpfr_mul (z_pre, z_pre, y, MPFR_RNDN);
+
           if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, prec1 - add, precz,
                                            rnd_mode)))
             break;
--- a/tests/tcan_round.c	2015-06-19 21:55:10.000000000 +0200
+++ b/tests/tcan_round.c	2016-02-16 14:55:39.863392138 +0100
@@ -1,4 +1,4 @@
-/* Test file for mpfr_can_round.
+/* Test file for mpfr_can_round and mpfr_round_p.
 
 Copyright 1999, 2001-2015 Free Software Foundation, Inc.
 Contributed by the AriC and Caramel projects, INRIA.
@@ -41,6 +41,8 @@
       /* avoid mpn_random which leaks memory */
       for (i = 0; i < n; i++)
         buf[i] = randlimb ();
+      /* force the number to be normalized */
+      buf[n - 1] |= MPFR_LIMB_HIGHBIT;
       p = randlimb() % ((n-1) * GMP_NUMB_BITS) + MPFR_PREC_MIN;
       err = p + randlimb () % GMP_NUMB_BITS;
       r1 = mpfr_round_p (buf, n, err, p);
@@ -57,11 +59,72 @@
     }
 }
 
+/* check x=2^i with precision px, error at most 1, and target precision prec */
+static void
+test_pow2 (mpfr_exp_t i, mpfr_prec_t px, mpfr_rnd_t r1, mpfr_rnd_t r2,
+           mpfr_prec_t prec)
+{
+  mpfr_t x;
+  int b, expected_b, b2;
+
+  mpfr_init2 (x, px);
+  mpfr_set_ui_2exp (x, 1, i, MPFR_RNDN);
+  b = !!mpfr_can_round (x, i+1, r1, r2, prec);
+  /* Note: If mpfr_can_round succeeds for both
+     (r1,r2) = (MPFR_RNDD,MPFR_RNDN) and
+     (r1,r2) = (MPFR_RNDU,MPFR_RNDN), then it should succeed for
+     (r1,r2) = (MPFR_RNDN,MPFR_RNDN). So, the condition on prec below
+     for r1 = MPFR_RNDN should be the most restrictive between those
+     for r1 = any directed rounding mode.
+     For r1 like MPFR_RNDA, the unrounded, unknown number may be anyone
+     in [2^i-1,i]. As both 2^i-1 and 2^i fit on i bits, one cannot round
+     in any precision >= i bits, hence the condition prec < i; prec = i-1
+     will work here for r2 = MPFR_RNDN thanks to the even-rounding rule
+     (and also with rounding ties away from zero). */
+  expected_b =
+    MPFR_IS_LIKE_RNDD (r1, MPFR_SIGN_POS) ?
+    (MPFR_IS_LIKE_RNDU (r2, MPFR_SIGN_POS) ? 0 : prec <= i) :
+    MPFR_IS_LIKE_RNDU (r1, MPFR_SIGN_POS) ?
+    (MPFR_IS_LIKE_RNDD (r2, MPFR_SIGN_POS) ? 0 : prec < i) :
+    (r2 != MPFR_RNDN ? 0 : prec < i);
+  /* We only require mpfr_can_round to return 1 when we can really
+     round, it is allowed to return 0 in some rare boundary cases,
+     for example when x = 2^k and the error is 0.25 ulp.
+     Note: if this changes in the future, the test could be improved by
+     removing the "&& expected_b == 0" below. */
+  if (b != expected_b && expected_b == 0)
+    {
+      printf ("Error for x=2^%d, px=%lu, err=%d, r1=%s, r2=%s, prec=%d\n",
+              (int) i, (unsigned long) px, (int) i + 1,
+              mpfr_print_rnd_mode (r1), mpfr_print_rnd_mode (r2), (int) prec);
+      printf ("Expected %d, got %d\n", expected_b, b);
+      exit (1);
+    }
+
+  if (r1 == MPFR_RNDN && r2 == MPFR_RNDZ)
+    {
+      /* Similar test to the one done in src/round_p.c
+         for MPFR_WANT_ASSERT >= 2. */
+      b2 = !!mpfr_round_p (MPFR_MANT(x), MPFR_LIMB_SIZE(x), i+1, prec);
+      if (b2 != b)
+        {
+          printf ("Error for x=2^%d, px=%lu, err=%d, prec=%d\n",
+                  (int) i, (unsigned long) px, (int) i + 1, (int) prec);
+          printf ("mpfr_can_round gave %d, mpfr_round_p gave %d\n", b, b2);
+          exit (1);
+        }
+    }
+
+  mpfr_clear (x);
+}
+
 int
 main (void)
 {
   mpfr_t x;
-  mpfr_prec_t i, j;
+  mpfr_prec_t i, j, k;
+  int r1, r2;
+  int n;
 
   tests_start_mpfr ();
 
@@ -111,12 +174,30 @@
   mpfr_set_str (x, "0.ff4ca619c76ba69", 16, MPFR_RNDZ);
   for (i = 30; i < 99; i++)
     for (j = 30; j < 99; j++)
-      {
-        int r1, r2;
-        for (r1 = 0; r1 < MPFR_RND_MAX ; r1++)
-          for (r2 = 0; r2 < MPFR_RND_MAX ; r2++)
-            mpfr_can_round (x, i, (mpfr_rnd_t) r1, (mpfr_rnd_t) r2, j); /* test for assertions */
-      }
+      for (r1 = 0; r1 < MPFR_RND_MAX; r1++)
+        for (r2 = 0; r2 < MPFR_RND_MAX; r2++)
+          {
+            /* test for assertions */
+            mpfr_can_round (x, i, (mpfr_rnd_t) r1, (mpfr_rnd_t) r2, j);
+          }
+
+  test_pow2 (32, 32, MPFR_RNDN, MPFR_RNDN, 32);
+  test_pow2 (174, 174, MPFR_RNDN, MPFR_RNDN, 174);
+  test_pow2 (174, 174, MPFR_RNDU, MPFR_RNDN, 174);
+  test_pow2 (176, 129, MPFR_RNDU, MPFR_RNDU, 174);
+  test_pow2 (176, 2, MPFR_RNDZ, MPFR_RNDZ, 174);
+  test_pow2 (176, 2, MPFR_RNDU, MPFR_RNDU, 176);
+
+  /* Tests for x = 2^i (E(x) = i+1) with error at most 1 = 2^0. */
+  for (n = 0; n < 100; n++)
+    {
+      i = (randlimb() % 200) + 4;
+      for (j = i - 2; j < i + 2; j++)
+        for (r1 = 0; r1 < MPFR_RND_MAX; r1++)
+          for (r2 = 0; r2 < MPFR_RND_MAX; r2++)
+            for (k = MPFR_PREC_MIN; k <= i + 2; k++)
+              test_pow2 (i, k, (mpfr_rnd_t) r1, (mpfr_rnd_t) r2, j);
+    }
 
   mpfr_clear (x);
 
--- a/tests/tdiv.c	2015-06-19 21:55:10.000000000 +0200
+++ b/tests/tdiv.c	2016-02-16 14:55:39.828391457 +0100
@@ -1099,6 +1099,69 @@
   mpfr_set_emax (old_emax);
 }
 
+/* Bug in mpfr_divhigh_n_basecase when all limbs of q (except the most
+   significant one) are B-1 where B=2^GMP_NUMB_BITS. Since we truncate
+   the divisor at each step, it might happen at some point that
+   (np[n-1],np[n-2]) > (d1,d0), and not only the equality.
+   Reported by Ricky Farr
+   <https://sympa.inria.fr/sympa/arc/mpfr/2015-10/msg00023.html>
+   To get a failure, a MPFR_DIVHIGH_TAB entry below the MPFR_DIV_THRESHOLD
+   limit must have a value 0. With most mparam.h files, this cannot occur. */
+static void
+test_20151023 (void)
+{
+  mpfr_prec_t p;
+  mpfr_t n, d, q, q0;
+  int inex, i;
+
+  for (p = GMP_NUMB_BITS; p <= 2000; p++)
+    {
+      mpfr_init2 (n, 2*p);
+      mpfr_init2 (d, p);
+      mpfr_init2 (q, p);
+      mpfr_init2 (q0, GMP_NUMB_BITS);
+
+      /* generate a random divisor of p bits */
+      mpfr_urandomb (d, RANDS);
+      /* generate a random quotient of GMP_NUMB_BITS bits */
+      mpfr_urandomb (q0, RANDS);
+      /* zero-pad the quotient to p bits */
+      inex = mpfr_prec_round (q0, p, MPFR_RNDN);
+      MPFR_ASSERTN(inex == 0);
+
+      for (i = 0; i < 3; i++)
+        {
+          /* i=0: try with the original quotient xxx000...000
+             i=1: try with the original quotient minus one ulp
+             i=2: try with the original quotient plus one ulp */
+          if (i == 1)
+            mpfr_nextbelow (q0);
+          else if (i == 2)
+            {
+              mpfr_nextabove (q0);
+              mpfr_nextabove (q0);
+            }
+
+          inex = mpfr_mul (n, d, q0, MPFR_RNDN);
+          MPFR_ASSERTN(inex == 0);
+          mpfr_nextabove (n);
+          mpfr_div (q, n, d, MPFR_RNDN);
+          MPFR_ASSERTN(mpfr_cmp (q, q0) == 0);
+
+          inex = mpfr_mul (n, d, q0, MPFR_RNDN);
+          MPFR_ASSERTN(inex == 0);
+          mpfr_nextbelow (n);
+          mpfr_div (q, n, d, MPFR_RNDN);
+          MPFR_ASSERTN(mpfr_cmp (q, q0) == 0);
+        }
+
+      mpfr_clear (n);
+      mpfr_clear (d);
+      mpfr_clear (q);
+      mpfr_clear (q0);
+    }
+}
+
 #define TEST_FUNCTION test_div
 #define TWO_ARGS
 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS)
@@ -1219,6 +1282,7 @@
   consistency ();
   test_20070603 ();
   test_20070628 ();
+  test_20151023 ();
   test_generic (2, 800, 50);
   test_extreme ();
 
--- a/tests/tfits.c	2015-06-19 21:55:10.000000000 +0200
+++ b/tests/tfits.c	2016-02-16 14:55:39.870392274 +0100
@@ -33,258 +33,225 @@
 #include "mpfr-intmax.h"
 #include "mpfr-test.h"
 
-#define ERROR1(N)                                               \
+#define FTEST_AUX(N,NOT,FCT)                                    \
   do                                                            \
     {                                                           \
-      printf("Error %d for rnd = %s and x = ", N,               \
-             mpfr_print_rnd_mode ((mpfr_rnd_t) r));             \
-      mpfr_dump(x);                                             \
-      exit(1);                                                  \
+      __gmpfr_flags = ex_flags;                                 \
+      if (NOT FCT (x, (mpfr_rnd_t) r))                          \
+        {                                                       \
+          printf ("Error %d for %s, rnd = %s and x = ",         \
+                  N, #FCT,                                      \
+                  mpfr_print_rnd_mode ((mpfr_rnd_t) r));        \
+          mpfr_dump (x);                                        \
+          exit (1);                                             \
+        }                                                       \
+      if (__gmpfr_flags != ex_flags)                            \
+        {                                                       \
+          unsigned int flags = __gmpfr_flags;                   \
+          printf ("Flags error %d for %s, rnd = %s and x = ",   \
+                  N, #FCT,                                      \
+                  mpfr_print_rnd_mode ((mpfr_rnd_t) r));        \
+          mpfr_dump(x);                                         \
+          printf ("Expected flags:");                           \
+          flags_out (ex_flags);                                 \
+          printf ("Got flags:     ");                           \
+          flags_out (flags);                                    \
+          exit (1);                                             \
+        }                                                       \
     }                                                           \
   while (0)
 
-static void check_intmax (void);
+#define FTEST(N,NOT,FCT)                                        \
+  do                                                            \
+    {                                                           \
+      mpfr_exp_t e;                                             \
+      FTEST_AUX (N,NOT,FCT);                                    \
+      if (MPFR_IS_SINGULAR (x))                                 \
+        break;                                                  \
+      e = mpfr_get_exp (x);                                     \
+      set_emin (e);                                             \
+      set_emax (e);                                             \
+      FTEST_AUX (N,NOT,FCT);                                    \
+      set_emin (emin);                                          \
+      set_emax (emax);                                          \
+    }                                                           \
+  while (0)
+
+#define CHECK_ALL(N,NOT)                                        \
+  do                                                            \
+    {                                                           \
+      FTEST (N, NOT, mpfr_fits_ulong_p);                        \
+      FTEST (N, NOT, mpfr_fits_slong_p);                        \
+      FTEST (N, NOT, mpfr_fits_uint_p);                         \
+      FTEST (N, NOT, mpfr_fits_sint_p);                         \
+      FTEST (N, NOT, mpfr_fits_ushort_p);                       \
+      FTEST (N, NOT, mpfr_fits_sshort_p);                       \
+    }                                                           \
+  while (0)
+
+#define CHECK_MAX(N,NOT)                                        \
+  do                                                            \
+    {                                                           \
+      FTEST (N, NOT, mpfr_fits_uintmax_p);                      \
+      FTEST (N, NOT, mpfr_fits_intmax_p);                       \
+    }                                                           \
+  while (0)
+
+/* V is a non-zero limit for the type (*_MIN for a signed type or *_MAX).
+ * If V is positive, then test V, V + 1/4, V + 3/4 and V + 1.
+ * If V is negative, then test V, V - 1/4, V - 3/4 and V - 1.
+ */
+#define CHECK_LIM(N,V,SET,FCT)                                  \
+  do                                                            \
+    {                                                           \
+      SET (x, V, MPFR_RNDN);                                    \
+      FTEST (N, !, FCT);                                        \
+      mpfr_set_si_2exp (y, (V) < 0 ? -1 : 1, -2, MPFR_RNDN);    \
+      mpfr_add (x, x, y, MPFR_RNDN);                            \
+      FTEST (N+1, (r == MPFR_RNDN ||                            \
+                   MPFR_IS_LIKE_RNDZ (r, (V) < 0)) ^ !!, FCT);  \
+      mpfr_add (x, x, y, MPFR_RNDN);                            \
+      mpfr_add (x, x, y, MPFR_RNDN);                            \
+      FTEST (N+3, MPFR_IS_LIKE_RNDZ (r, (V) < 0) ^ !!, FCT);    \
+      mpfr_add (x, x, y, MPFR_RNDN);                            \
+      FTEST (N+4, !!, FCT);                                     \
+    }                                                           \
+  while (0)
 
 int
 main (void)
 {
+  mpfr_exp_t emin, emax;
   mpfr_t x, y;
-  int i, r;
+  unsigned int flags[2] = { 0, MPFR_FLAGS_ALL }, ex_flags;
+  int i, r, fi;
 
   tests_start_mpfr ();
 
-  mpfr_init2 (x, 256);
+  emin = mpfr_get_emin ();
+  emax = mpfr_get_emax ();
+
+  mpfr_init2 (x, sizeof (unsigned long) * CHAR_BIT + 2);
   mpfr_init2 (y, 8);
 
   RND_LOOP (r)
-    {
-
-      /* Check NAN */
-      mpfr_set_nan (x);
-      if (mpfr_fits_ulong_p (x, (mpfr_rnd_t) r))
-        ERROR1 (1);
-      if (mpfr_fits_slong_p (x, (mpfr_rnd_t) r))
-        ERROR1 (2);
-      if (mpfr_fits_uint_p (x, (mpfr_rnd_t) r))
-        ERROR1 (3);
-      if (mpfr_fits_sint_p (x, (mpfr_rnd_t) r))
-        ERROR1 (4);
-      if (mpfr_fits_ushort_p (x, (mpfr_rnd_t) r))
-        ERROR1 (5);
-      if (mpfr_fits_sshort_p (x, (mpfr_rnd_t) r))
-        ERROR1 (6);
-
-      /* Check INF */
-      mpfr_set_inf (x, 1);
-      if (mpfr_fits_ulong_p (x, (mpfr_rnd_t) r))
-        ERROR1 (7);
-      if (mpfr_fits_slong_p (x, (mpfr_rnd_t) r))
-        ERROR1 (8);
-      if (mpfr_fits_uint_p (x, (mpfr_rnd_t) r))
-        ERROR1 (9);
-      if (mpfr_fits_sint_p (x, (mpfr_rnd_t) r))
-        ERROR1 (10);
-      if (mpfr_fits_ushort_p (x, (mpfr_rnd_t) r))
-        ERROR1 (11);
-      if (mpfr_fits_sshort_p (x, (mpfr_rnd_t) r))
-        ERROR1 (12);
-
-      /* Check Zero */
-      MPFR_SET_ZERO (x);
-      if (!mpfr_fits_ulong_p (x, (mpfr_rnd_t) r))
-        ERROR1 (13);
-      if (!mpfr_fits_slong_p (x, (mpfr_rnd_t) r))
-        ERROR1 (14);
-      if (!mpfr_fits_uint_p (x, (mpfr_rnd_t) r))
-        ERROR1 (15);
-      if (!mpfr_fits_sint_p (x, (mpfr_rnd_t) r))
-        ERROR1 (16);
-      if (!mpfr_fits_ushort_p (x, (mpfr_rnd_t) r))
-        ERROR1 (17);
-      if (!mpfr_fits_sshort_p (x, (mpfr_rnd_t) r))
-        ERROR1 (18);
-
-      /* Check small positive op */
-      mpfr_set_str1 (x, "1@-1");
-      if (!mpfr_fits_ulong_p (x, (mpfr_rnd_t) r))
-        ERROR1 (19);
-      if (!mpfr_fits_slong_p (x, (mpfr_rnd_t) r))
-        ERROR1 (20);
-      if (!mpfr_fits_uint_p (x, (mpfr_rnd_t) r))
-        ERROR1 (21);
-      if (!mpfr_fits_sint_p (x, (mpfr_rnd_t) r))
-        ERROR1 (22);
-      if (!mpfr_fits_ushort_p (x, (mpfr_rnd_t) r))
-        ERROR1 (23);
-      if (!mpfr_fits_sshort_p (x, (mpfr_rnd_t) r))
-        ERROR1 (24);
-
-      /* Check 17 */
-      mpfr_set_ui (x, 17, MPFR_RNDN);
-      if (!mpfr_fits_ulong_p (x, (mpfr_rnd_t) r))
-        ERROR1 (25);
-      if (!mpfr_fits_slong_p (x, (mpfr_rnd_t) r))
-        ERROR1 (26);
-      if (!mpfr_fits_uint_p (x, (mpfr_rnd_t) r))
-        ERROR1 (27);
-      if (!mpfr_fits_sint_p (x, (mpfr_rnd_t) r))
-        ERROR1 (28);
-      if (!mpfr_fits_ushort_p (x, (mpfr_rnd_t) r))
-        ERROR1 (29);
-      if (!mpfr_fits_sshort_p (x, (mpfr_rnd_t) r))
-        ERROR1 (30);
-
-      /* Check all other values */
-      mpfr_set_ui (x, ULONG_MAX, MPFR_RNDN);
-      mpfr_mul_2exp (x, x, 1, MPFR_RNDN);
-      if (mpfr_fits_ulong_p (x, (mpfr_rnd_t) r))
-        ERROR1 (31);
-      if (mpfr_fits_slong_p (x, (mpfr_rnd_t) r))
-        ERROR1 (32);
-      mpfr_mul_2exp (x, x, 40, MPFR_RNDN);
-      if (mpfr_fits_ulong_p (x, (mpfr_rnd_t) r))
-        ERROR1 (33);
-      if (mpfr_fits_uint_p (x, (mpfr_rnd_t) r))
-        ERROR1 (34);
-      if (mpfr_fits_sint_p (x, (mpfr_rnd_t) r))
-        ERROR1 (35);
-      if (mpfr_fits_ushort_p (x, (mpfr_rnd_t) r))
-        ERROR1 (36);
-      if (mpfr_fits_sshort_p (x, (mpfr_rnd_t) r))
-        ERROR1 (37);
-
-      mpfr_set_ui (x, ULONG_MAX, MPFR_RNDN);
-      if (!mpfr_fits_ulong_p (x, (mpfr_rnd_t) r))
-        ERROR1 (38);
-      mpfr_set_ui (x, LONG_MAX, MPFR_RNDN);
-      if (!mpfr_fits_slong_p (x, (mpfr_rnd_t) r))
-        ERROR1 (39);
-      mpfr_set_ui (x, UINT_MAX, MPFR_RNDN);
-      if (!mpfr_fits_uint_p (x, (mpfr_rnd_t) r))
-        ERROR1 (40);
-      mpfr_set_ui (x, INT_MAX, MPFR_RNDN);
-      if (!mpfr_fits_sint_p (x, (mpfr_rnd_t) r))
-        ERROR1 (41);
-      mpfr_set_ui (x, USHRT_MAX, MPFR_RNDN);
-      if (!mpfr_fits_ushort_p (x, (mpfr_rnd_t) r))
-        ERROR1 (42);
-      mpfr_set_ui (x, SHRT_MAX, MPFR_RNDN);
-      if (!mpfr_fits_sshort_p (x, (mpfr_rnd_t) r))
-        ERROR1 (43);
-
-      mpfr_set_si (x, 1, MPFR_RNDN);
-      if (!mpfr_fits_sint_p (x, (mpfr_rnd_t) r))
-        ERROR1 (44);
-      if (!mpfr_fits_sshort_p (x, (mpfr_rnd_t) r))
-        ERROR1 (45);
-
-      /* Check negative op */
-      for (i = 1; i <= 4; i++)
-        {
-          int inv;
-
-          mpfr_set_si_2exp (x, -i, -2, MPFR_RNDN);
-          mpfr_rint (y, x, (mpfr_rnd_t) r);
-          inv = MPFR_NOTZERO (y);
-          if (!mpfr_fits_ulong_p (x, (mpfr_rnd_t) r) ^ inv)
-            ERROR1 (46);
-          if (!mpfr_fits_slong_p (x, (mpfr_rnd_t) r))
-            ERROR1 (47);
-          if (!mpfr_fits_uint_p (x, (mpfr_rnd_t) r) ^ inv)
-            ERROR1 (48);
-          if (!mpfr_fits_sint_p (x, (mpfr_rnd_t) r))
-            ERROR1 (49);
-          if (!mpfr_fits_ushort_p (x, (mpfr_rnd_t) r) ^ inv)
-            ERROR1 (50);
-          if (!mpfr_fits_sshort_p (x, (mpfr_rnd_t) r))
-            ERROR1 (51);
-        }
-    }
-
-  mpfr_clear (x);
-  mpfr_clear (y);
+    for (fi = 0; fi < numberof (flags); fi++)
+      {
+        ex_flags = flags[fi];
+
+        /* Check NaN */
+        mpfr_set_nan (x);
+        CHECK_ALL (1, !!);
+
+        /* Check +Inf */
+        mpfr_set_inf (x, 1);
+        CHECK_ALL (2, !!);
+
+        /* Check -Inf */
+        mpfr_set_inf (x, -1);
+        CHECK_ALL (3, !!);
+
+        /* Check +0 */
+        mpfr_set_zero (x, 1);
+        CHECK_ALL (4, !);
+
+        /* Check -0 */
+        mpfr_set_zero (x, -1);
+        CHECK_ALL (5, !);
+
+        /* Check small positive op */
+        mpfr_set_str1 (x, "1@-1");
+        CHECK_ALL (6, !);
+
+        /* Check 17 */
+        mpfr_set_ui (x, 17, MPFR_RNDN);
+        CHECK_ALL (7, !);
+
+        /* Check large values (no fit) */
+        mpfr_set_ui (x, ULONG_MAX, MPFR_RNDN);
+        mpfr_mul_2exp (x, x, 1, MPFR_RNDN);
+        CHECK_ALL (8, !!);
+        mpfr_mul_2exp (x, x, 40, MPFR_RNDN);
+        CHECK_ALL (9, !!);
+
+        /* Check a non-integer number just below a power of two. */
+        mpfr_set_ui_2exp (x, 255, -2, MPFR_RNDN);
+        CHECK_ALL (10, !);
+
+        /* Check the limits of the types (except 0 for unsigned types) */
+        CHECK_LIM (20, ULONG_MAX, mpfr_set_ui, mpfr_fits_ulong_p);
+        CHECK_LIM (30, LONG_MAX, mpfr_set_si, mpfr_fits_slong_p);
+        CHECK_LIM (35, LONG_MIN, mpfr_set_si, mpfr_fits_slong_p);
+        CHECK_LIM (40, UINT_MAX, mpfr_set_ui, mpfr_fits_uint_p);
+        CHECK_LIM (50, INT_MAX, mpfr_set_si, mpfr_fits_sint_p);
+        CHECK_LIM (55, INT_MIN, mpfr_set_si, mpfr_fits_sint_p);
+        CHECK_LIM (60, USHRT_MAX, mpfr_set_ui, mpfr_fits_ushort_p);
+        CHECK_LIM (70, SHRT_MAX, mpfr_set_si, mpfr_fits_sshort_p);
+        CHECK_LIM (75, SHRT_MIN, mpfr_set_si, mpfr_fits_sshort_p);
+
+        /* Check negative op */
+        for (i = 1; i <= 4; i++)
+          {
+            int inv;
+
+            mpfr_set_si_2exp (x, -i, -2, MPFR_RNDN);
+            mpfr_rint (y, x, (mpfr_rnd_t) r);
+            inv = MPFR_NOTZERO (y);
+            FTEST (80, inv ^ !, mpfr_fits_ulong_p);
+            FTEST (81,       !, mpfr_fits_slong_p);
+            FTEST (82, inv ^ !, mpfr_fits_uint_p);
+            FTEST (83,       !, mpfr_fits_sint_p);
+            FTEST (84, inv ^ !, mpfr_fits_ushort_p);
+            FTEST (85,       !, mpfr_fits_sshort_p);
+          }
+      }
 
-  check_intmax ();
-
-  tests_end_mpfr ();
-  return 0;
-}
-
-static void
-check_intmax (void)
-{
 #ifdef _MPFR_H_HAVE_INTMAX_T
-  mpfr_t x, y;
-  int i, r;
 
-  mpfr_init2 (x, sizeof (uintmax_t) * CHAR_BIT);
-  mpfr_init2 (y, 8);
+  mpfr_set_prec (x, sizeof (uintmax_t) * CHAR_BIT + 2);
 
   RND_LOOP (r)
     {
-      /* Check NAN */
+      /* Check NaN */
       mpfr_set_nan (x);
-      if (mpfr_fits_uintmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (52);
-      if (mpfr_fits_intmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (53);
+      CHECK_MAX (1, !!);
 
-      /* Check INF */
+      /* Check +Inf */
       mpfr_set_inf (x, 1);
-      if (mpfr_fits_uintmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (54);
-      if (mpfr_fits_intmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (55);
-
-      /* Check Zero */
-      MPFR_SET_ZERO (x);
-      if (!mpfr_fits_uintmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (56);
-      if (!mpfr_fits_intmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (57);
+      CHECK_MAX (2, !!);
 
-      /* Check positive small op */
+      /* Check -Inf */
+      mpfr_set_inf (x, -1);
+      CHECK_MAX (3, !!);
+
+      /* Check +0 */
+      mpfr_set_zero (x, 1);
+      CHECK_MAX (4, !);
+
+      /* Check -0 */
+      mpfr_set_zero (x, -1);
+      CHECK_MAX (5, !);
+
+      /* Check small positive op */
       mpfr_set_str1 (x, "1@-1");
-      if (!mpfr_fits_uintmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (58);
-      if (!mpfr_fits_intmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (59);
+      CHECK_MAX (6, !);
 
       /* Check 17 */
       mpfr_set_ui (x, 17, MPFR_RNDN);
-      if (!mpfr_fits_uintmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (60);
-      if (!mpfr_fits_intmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (61);
+      CHECK_MAX (7, !);
 
       /* Check hugest */
       mpfr_set_ui_2exp (x, 42, sizeof (uintmax_t) * 32, MPFR_RNDN);
-      if (mpfr_fits_uintmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (62);
-      if (mpfr_fits_intmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (63);
-
-      /* Check all other values */
-      mpfr_set_uj (x, MPFR_UINTMAX_MAX, MPFR_RNDN);
-      mpfr_add_ui (x, x, 1, MPFR_RNDN);
-      if (mpfr_fits_uintmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (64);
-      mpfr_set_uj (x, MPFR_UINTMAX_MAX, MPFR_RNDN);
-      if (!mpfr_fits_uintmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (65);
-      mpfr_set_sj (x, MPFR_INTMAX_MAX, MPFR_RNDN);
-      mpfr_add_ui (x, x, 1, MPFR_RNDN);
-      if (mpfr_fits_intmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (66);
-      mpfr_set_sj (x, MPFR_INTMAX_MAX, MPFR_RNDN);
-      if (!mpfr_fits_intmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (67);
-      mpfr_set_sj (x, MPFR_INTMAX_MIN, MPFR_RNDN);
-      if (!mpfr_fits_intmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (68);
-      mpfr_sub_ui (x, x, 1, MPFR_RNDN);
-      if (mpfr_fits_intmax_p (x, (mpfr_rnd_t) r))
-        ERROR1 (69);
+      CHECK_MAX (8, !!);
+
+      /* Check a non-integer number just below a power of two. */
+      mpfr_set_ui_2exp (x, 255, -2, MPFR_RNDN);
+      CHECK_MAX (10, !);
+
+      /* Check the limits of the types (except 0 for uintmax_t) */
+      CHECK_LIM (20, MPFR_UINTMAX_MAX, mpfr_set_uj, mpfr_fits_uintmax_p);
+      CHECK_LIM (30, MPFR_INTMAX_MAX, mpfr_set_sj, mpfr_fits_intmax_p);
+      CHECK_LIM (35, MPFR_INTMAX_MIN, mpfr_set_sj, mpfr_fits_intmax_p);
 
       /* Check negative op */
       for (i = 1; i <= 4; i++)
@@ -294,14 +261,16 @@
           mpfr_set_si_2exp (x, -i, -2, MPFR_RNDN);
           mpfr_rint (y, x, (mpfr_rnd_t) r);
           inv = MPFR_NOTZERO (y);
-          if (!mpfr_fits_uintmax_p (x, (mpfr_rnd_t) r) ^ inv)
-            ERROR1 (70);
-          if (!mpfr_fits_intmax_p (x, (mpfr_rnd_t) r))
-            ERROR1 (71);
+          FTEST (80, inv ^ !, mpfr_fits_uintmax_p);
+          FTEST (81,       !, mpfr_fits_intmax_p);
         }
     }
 
+#endif  /* _MPFR_H_HAVE_INTMAX_T */
+
   mpfr_clear (x);
   mpfr_clear (y);
-#endif
+
+  tests_end_mpfr ();
+  return 0;
 }
--- a/tests/tfrexp.c	2015-06-19 21:55:10.000000000 +0200
+++ b/tests/tfrexp.c	2016-02-16 14:55:39.824391379 +0100
@@ -129,12 +129,115 @@
   mpfr_clear (x);
 }
 
+static void check1 (void)
+{
+  mpfr_exp_t emin, emax, e;
+  mpfr_t x, y1, y2;
+  int r, neg, red;
+
+  emin = mpfr_get_emin ();
+  emax = mpfr_get_emax ();
+  set_emin (MPFR_EMIN_MIN);
+  set_emax (MPFR_EMAX_MAX);
+
+  mpfr_init2 (x, 7);
+  mpfr_inits2 (4, y1, y2, (mpfr_ptr) 0);
+
+  mpfr_set_ui_2exp (x, 1, -2, MPFR_RNDN);
+  while (mpfr_regular_p (x))
+    {
+      /* Test the exponents up to 3 and with the maximum exponent
+         (to check potential intermediate overflow). */
+      if (MPFR_GET_EXP (x) == 4)
+        mpfr_set_exp (x, MPFR_EMAX_MAX);
+      e = MPFR_GET_EXP (x);
+      for (neg = 0; neg < 2; neg++)
+        {
+          RND_LOOP (r)
+            {
+              int inex1, inex2;
+              mpfr_exp_t e1, e2;
+              unsigned int flags1, flags2;
+
+              for (red = 0; red < 2; red++)
+                {
+                  if (red)
+                    {
+                      /* e1: exponent of the rounded value of x. */
+                      MPFR_ASSERTN (e1 == e || e1 == e + 1);
+                      set_emin (e);
+                      set_emax (e);
+                      mpfr_clear_flags ();
+                      inex1 = e1 < 0 ?
+                        mpfr_mul_2ui (y1, x, -e1, (mpfr_rnd_t) r) :
+                        mpfr_div_2ui (y1, x, e1, (mpfr_rnd_t) r);
+                      flags1 = __gmpfr_flags;
+                    }
+                  else
+                    {
+                      inex1 = mpfr_set (y1, x, (mpfr_rnd_t) r);
+                      e1 = MPFR_IS_INF (y1) ? e + 1 : MPFR_GET_EXP (y1);
+                      flags1 = inex1 != 0 ? MPFR_FLAGS_INEXACT : 0;
+                    }
+                  mpfr_clear_flags ();
+                  inex2 = mpfr_frexp (&e2, y2, x, (mpfr_rnd_t) r);
+                  flags2 = __gmpfr_flags;
+                  set_emin (MPFR_EMIN_MIN);
+                  set_emax (MPFR_EMAX_MAX);
+                  if ((!red || e == 0) &&
+                      (! mpfr_regular_p (y2) || MPFR_GET_EXP (y2) != 0))
+                    {
+                      printf ("Error in check1 for %s, red = %d, x = ",
+                              mpfr_print_rnd_mode ((mpfr_rnd_t) r), red);
+                      mpfr_dump (x);
+                      printf ("Expected 1/2 <= |y| < 1, got y = ");
+                      mpfr_dump (y2);
+                      exit (1);
+                    }
+                  if (!red)
+                    {
+                      if (e2 > 0)
+                        mpfr_mul_2ui (y2, y2, e2, MPFR_RNDN);
+                      else if (e2 < 0)
+                        mpfr_div_2ui (y2, y2, -e2, MPFR_RNDN);
+                    }
+                  if (! (SAME_SIGN (inex1, inex2) &&
+                         mpfr_equal_p (y1, y2) &&
+                         flags1 == flags2))
+                    {
+                      printf ("Error in check1 for %s, red = %d, x = ",
+                              mpfr_print_rnd_mode ((mpfr_rnd_t) r), red);
+                      mpfr_dump (x);
+                      printf ("Expected y1 = ");
+                      mpfr_dump (y1);
+                      printf ("Got      y2 = ");
+                      mpfr_dump (y2);
+                      printf ("Expected inex ~= %d, got %d\n", inex1, inex2);
+                      printf ("Expected flags:");
+                      flags_out (flags1);
+                      printf ("Got flags:     ");
+                      flags_out (flags2);
+                      exit (1);
+                    }
+                }
+            }
+          mpfr_neg (x, x, MPFR_RNDN);
+        }
+      mpfr_nextabove (x);
+    }
+
+  mpfr_clears (x, y1, y2, (mpfr_ptr) 0);
+  set_emin (emin);
+  set_emax (emax);
+}
+
 int
 main (int argc, char *argv[])
 {
   tests_start_mpfr ();
 
   check_special ();
+  check1 ();
 
   tests_end_mpfr ();
   return 0;
--- a/tests/tj0.c	2015-06-19 21:55:10.000000000 +0200
+++ b/tests/tj0.c	2016-02-16 14:55:39.831391515 +0100
@@ -99,6 +99,18 @@
   mpfr_j0 (y, x, MPFR_RNDN);
   MPFR_ASSERTN (! mpfr_nan_p (y) && mpfr_cmp_ui_2exp (y, 41, -11) == 0);
 
+  /* Bug reported by Fredrik Johansson on 19 Jan 2016 */
+  mpfr_set_prec (x, 53);
+  mpfr_set_str (x, "0x4.3328p+0", 0, MPFR_RNDN);
+  mpfr_set_prec (y, 2);
+  mpfr_j0 (y, x, MPFR_RNDD);
+  /* y should be -0.5 */
+  MPFR_ASSERTN (! mpfr_nan_p (y) && mpfr_cmp_si_2exp (y, -1, -1) == 0);
+  mpfr_set_prec (y, 3);
+  mpfr_j0 (y, x, MPFR_RNDD);
+  /* y should be -0.4375 */
+  MPFR_ASSERTN (! mpfr_nan_p (y) && mpfr_cmp_si_2exp (y, -7, -4) == 0);
+
   /* Case for which s = 0 in mpfr_jn */
   mpfr_set_prec (x, 44);
   mpfr_set_prec (y, 44);
--- a/tests/tlngamma.c	2015-06-19 21:55:10.000000000 +0200
+++ b/tests/tlngamma.c	2016-02-16 14:55:39.809391087 +0100
@@ -33,7 +33,7 @@
 special (void)
 {
   mpfr_t x, y;
-  int inex;
+  int i, inex;
 
   mpfr_init (x);
   mpfr_init (y);
@@ -46,25 +46,29 @@
       exit (1);
     }
 
-  mpfr_set_inf (x, -1);
+  mpfr_set_inf (x, 1);
+  mpfr_clear_flags ();
   mpfr_lngamma (y, x, MPFR_RNDN);
-  if (!mpfr_nan_p (y))
+  if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 || __gmpfr_flags != 0)
     {
-      printf ("Error for lngamma(-Inf)\n");
+      printf ("Error for lngamma(+Inf)\n");
       exit (1);
     }
 
-  mpfr_set_inf (x, 1);
+  mpfr_set_inf (x, -1);
+  mpfr_clear_flags ();
   mpfr_lngamma (y, x, MPFR_RNDN);
-  if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
+  if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 || __gmpfr_flags != 0)
     {
-      printf ("Error for lngamma(+Inf)\n");
+      printf ("Error for lngamma(-Inf)\n");
       exit (1);
     }
 
   mpfr_set_ui (x, 0, MPFR_RNDN);
+  mpfr_clear_flags ();
   mpfr_lngamma (y, x, MPFR_RNDN);
-  if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
+  if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 ||
+      __gmpfr_flags != MPFR_FLAGS_DIVBY0)
     {
       printf ("Error for lngamma(+0)\n");
       exit (1);
@@ -72,32 +76,58 @@
 
   mpfr_set_ui (x, 0, MPFR_RNDN);
   mpfr_neg (x, x, MPFR_RNDN);
+  mpfr_clear_flags ();
   mpfr_lngamma (y, x, MPFR_RNDN);
-  if (!mpfr_nan_p (y))
+  if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 ||
+      __gmpfr_flags != MPFR_FLAGS_DIVBY0)
     {
       printf ("Error for lngamma(-0)\n");
       exit (1);
     }
 
   mpfr_set_ui (x, 1, MPFR_RNDN);
+  mpfr_clear_flags ();
   mpfr_lngamma (y, x, MPFR_RNDN);
-  if (MPFR_IS_NAN (y) || mpfr_cmp_ui (y, 0) || MPFR_IS_NEG (y))
+  if (mpfr_cmp_ui0 (y, 0) || MPFR_IS_NEG (y))
     {
       printf ("Error for lngamma(1)\n");
       exit (1);
     }
 
-  mpfr_set_si (x, -1, MPFR_RNDN);
-  mpfr_lngamma (y, x, MPFR_RNDN);
-  if (!mpfr_nan_p (y))
+  for (i = 1; i <= 5; i++)
     {
-      printf ("Error for lngamma(-1)\n");
-      exit (1);
+      int c;
+
+      mpfr_set_si (x, -i, MPFR_RNDN);
+      mpfr_clear_flags ();
+      mpfr_lngamma (y, x, MPFR_RNDN);
+      if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 ||
+          __gmpfr_flags != MPFR_FLAGS_DIVBY0)
+        {
+          printf ("Error for lngamma(-%d)\n", i);
+          exit (1);
+        }
+      if (i & 1)
+        {
+          mpfr_nextabove (x);
+          c = '+';
+        }
+      else
+        {
+          mpfr_nextbelow (x);
+          c = '-';
+        }
+      mpfr_lngamma (y, x, MPFR_RNDN);
+      if (!mpfr_nan_p (y))
+        {
+          printf ("Error for lngamma(-%d%cepsilon)\n", i, c);
+          exit (1);
+        }
     }
 
   mpfr_set_ui (x, 2, MPFR_RNDN);
   mpfr_lngamma (y, x, MPFR_RNDN);
-  if (MPFR_IS_NAN (y) || mpfr_cmp_ui (y, 0) || MPFR_IS_NEG (y))
+  if (mpfr_cmp_ui0 (y, 0) || MPFR_IS_NEG (y))
     {
       printf ("Error for lngamma(2)\n");
       exit (1);
@@ -127,7 +157,7 @@
   mpfr_set_str (x, CHECK_X2, 10, MPFR_RNDN);
   mpfr_lngamma (y, x, MPFR_RNDN);
   mpfr_set_str (x, CHECK_Y2, 10, MPFR_RNDN);
-  if (MPFR_IS_NAN (y) || mpfr_cmp (y, x))
+  if (mpfr_cmp0 (y, x))
     {
       printf ("mpfr_lngamma("CHECK_X2") is wrong:\n"
               "expected ");
@@ -143,7 +173,7 @@
   mpfr_lngamma (y, x, MPFR_RNDU);
   mpfr_set_prec (x, 175);
   mpfr_set_str_binary (x, "0.1010001100011101101011001101110010100001000001000001110011000001101100001111001001000101011011100100010101011110100111110101010100010011010010000101010111001100011000101111E7");
-  if (MPFR_IS_NAN (y) || mpfr_cmp (x, y))
+  if (mpfr_cmp0 (x, y))
     {
       printf ("Error in mpfr_lngamma (1)\n");
       exit (1);
@@ -155,7 +185,7 @@
   mpfr_lngamma (x, y, MPFR_RNDZ);
   mpfr_set_prec (y, 21);
   mpfr_set_str_binary (y, "0.111000101000001100101E9");
-  if (MPFR_IS_NAN (x) || mpfr_cmp (x, y))
+  if (mpfr_cmp0 (x, y))
     {
       printf ("Error in mpfr_lngamma (120)\n");
       printf ("Expected "); mpfr_print_binary (y); puts ("");
@@ -169,7 +199,7 @@
   inex = mpfr_lngamma (y, x, MPFR_RNDN);
   mpfr_set_prec (x, 206);
   mpfr_set_str_binary (x, "0.10000111011000000011100010101001100110001110000111100011000100100110110010001011011110101001111011110110000001010100111011010000000011100110110101100111000111010011110010000100010111101010001101000110101001E13");
-  if (MPFR_IS_NAN (y) || mpfr_cmp (x, y))
+  if (mpfr_cmp0 (x, y))
     {
       printf ("Error in mpfr_lngamma (768)\n");
       exit (1);
@@ -185,7 +215,7 @@
   mpfr_set_str_binary (x, "0.1100E-66");
   mpfr_lngamma (y, x, MPFR_RNDN);
   mpfr_set_str_binary (x, "0.1100E6");
-  if (MPFR_IS_NAN (y) || mpfr_cmp (x, y))
+  if (mpfr_cmp0 (x, y))
     {
       printf ("Error for lngamma(0.1100E-66)\n");
       exit (1);
@@ -199,7 +229,7 @@
   mpfr_lngamma (y, x, MPFR_RNDN);
   mpfr_set_prec (x, 32);
   mpfr_set_str_binary (x, "-0.10001000111011111011000010100010E207");
-  if (MPFR_IS_NAN (y) || mpfr_cmp (x, y))
+  if (mpfr_cmp0 (x, y))
     {
       printf ("Error for lngamma(-2^199+0.5)\n");
       printf ("Got        ");
--- a/tests/tmul_2exp.c	2015-06-19 21:55:10.000000000 +0200
+++ b/tests/tmul_2exp.c	2016-02-16 14:55:39.818391262 +0100
@@ -50,77 +50,82 @@
 {
   mpfr_t x, y, z1, z2;
   mpfr_exp_t emin;
-  int i, k;
+  int i, k, s;
   int prec;
   int rnd;
   int div;
   int inex1, inex2;
   unsigned int flags1, flags2;
 
-  /* Test mul_2si(x, e - k), div_2si(x, k - e) and div_2ui(x, k - e)
-   * with emin = e, x = 1 + i/16, i in { -1, 0, 1 }, and k = 1 to 4,
-   * by comparing the result with the one of a simple division.
+  /* Test mul_2si(x, e - k), div_2si(x, k - e) and div_2ui(x, k - e) with
+   * emin = e, x = s * (1 + i/16), i in { -1, 0, 1 }, s in { -1, 1 }, and
+   * k = 1 to 4, by comparing the result with the one of a simple division.
    */
   emin = mpfr_get_emin ();
   set_emin (e);
   mpfr_inits2 (8, x, y, (mpfr_ptr) 0);
   for (i = 15; i <= 17; i++)
-    {
-      inex1 = mpfr_set_ui_2exp (x, i, -4, MPFR_RNDN);
-      MPFR_ASSERTN (inex1 == 0);
-      for (prec = 6; prec >= 3; prec -= 3)
-        {
-          mpfr_inits2 (prec, z1, z2, (mpfr_ptr) 0);
-          RND_LOOP (rnd)
-            for (k = 1; k <= 4; k++)
-              {
-                /* The following one is assumed to be correct. */
-                inex1 = mpfr_mul_2si (y, x, e, MPFR_RNDN);
-                MPFR_ASSERTN (inex1 == 0);
-                inex1 = mpfr_set_ui (z1, 1 << k, MPFR_RNDN);
-                MPFR_ASSERTN (inex1 == 0);
-                mpfr_clear_flags ();
-                /* Do not use mpfr_div_ui to avoid the optimization
-                   by mpfr_div_2si. */
-                inex1 = mpfr_div (z1, y, z1, (mpfr_rnd_t) rnd);
-                flags1 = __gmpfr_flags;
-
-              for (div = 0; div <= 2; div++)
+    for (s = 1; s >= -1; s -= 2)
+      {
+        inex1 = mpfr_set_si_2exp (x, s * i, -4, MPFR_RNDN);
+        MPFR_ASSERTN (inex1 == 0);
+        for (prec = 6; prec >= 3; prec -= 3)
+          {
+            mpfr_inits2 (prec, z1, z2, (mpfr_ptr) 0);
+            RND_LOOP (rnd)
+              for (k = 1; k <= 4; k++)
                 {
+                  /* The following one is assumed to be correct. */
+                  inex1 = mpfr_mul_2si (y, x, e, MPFR_RNDN);
+                  MPFR_ASSERTN (inex1 == 0);
+                  inex1 = mpfr_set_ui (z1, 1 << k, MPFR_RNDN);
+                  MPFR_ASSERTN (inex1 == 0);
                   mpfr_clear_flags ();
-                  inex2 = div == 0 ?
-                    mpfr_mul_2si (z2, x, e - k, (mpfr_rnd_t) rnd) : div == 1 ?
-                    mpfr_div_2si (z2, x, k - e, (mpfr_rnd_t) rnd) :
-                    mpfr_div_2ui (z2, x, k - e, (mpfr_rnd_t) rnd);
-                  flags2 = __gmpfr_flags;
-                  if (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
-                      mpfr_equal_p (z1, z2))
-                    continue;
-                  printf ("Error in underflow(");
-                  if (e == MPFR_EMIN_MIN)
-                    printf ("MPFR_EMIN_MIN");
-                  else if (e == emin)
-                    printf ("default emin");
-                  else if (e >= LONG_MIN)
-                    printf ("%ld", (long) e);
-                  else
-                    printf ("<LONG_MIN");
-                  printf (") with mpfr_%s,\nx = %d/16, prec = %d, k = %d, "
-                          "%s\n", div == 0 ? "mul_2si" : div == 1 ?
-                          "div_2si" : "div_2ui", i, prec, k,
-                          mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
-                  printf ("Expected ");
-                  mpfr_out_str (stdout, 16, 0, z1, MPFR_RNDN);
-                  printf (", inex = %d, flags = %u\n", SIGN (inex1), flags1);
-                  printf ("Got      ");
-                  mpfr_out_str (stdout, 16, 0, z2, MPFR_RNDN);
-                  printf (", inex = %d, flags = %u\n", SIGN (inex2), flags2);
-                  exit (1);
-                }  /* div */
-              }  /* k */
-          mpfr_clears (z1, z2, (mpfr_ptr) 0);
-        }  /* prec */
-    }  /* i */
+                  /* Do not use mpfr_div_ui to avoid the optimization
+                     by mpfr_div_2si. */
+                  inex1 = mpfr_div (z1, y, z1, (mpfr_rnd_t) rnd);
+                  flags1 = __gmpfr_flags;
+
+                  for (div = 0; div <= 2; div++)
+                    {
+                      mpfr_clear_flags ();
+                      inex2 =
+                        div == 0 ?
+                        mpfr_mul_2si (z2, x, e - k, (mpfr_rnd_t) rnd) :
+                        div == 1 ?
+                        mpfr_div_2si (z2, x, k - e, (mpfr_rnd_t) rnd) :
+                        mpfr_div_2ui (z2, x, k - e, (mpfr_rnd_t) rnd);
+                      flags2 = __gmpfr_flags;
+                      if (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
+                          mpfr_equal_p (z1, z2))
+                        continue;
+                      printf ("Error in underflow(");
+                      if (e == MPFR_EMIN_MIN)
+                        printf ("MPFR_EMIN_MIN");
+                      else if (e == emin)
+                        printf ("default emin");
+                      else if (e >= LONG_MIN)
+                        printf ("%ld", (long) e);
+                      else
+                        printf ("<LONG_MIN");
+                      printf (") with mpfr_%s,\nx = %d/16, prec = %d, k = %d,"
+                              " %s\n", div == 0 ? "mul_2si" : div == 1 ?
+                              "div_2si" : "div_2ui", s * i, prec, k,
+                              mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
+                      printf ("Expected ");
+                      mpfr_out_str (stdout, 16, 0, z1, MPFR_RNDN);
+                      printf (", inex = %d, flags = %u\n",
+                              SIGN (inex1), flags1);
+                      printf ("Got      ");
+                      mpfr_out_str (stdout, 16, 0, z2, MPFR_RNDN);
+                      printf (", inex = %d, flags = %u\n",
+                              SIGN (inex2), flags2);
+                      exit (1);
+                    }  /* div */
+                }  /* k */
+            mpfr_clears (z1, z2, (mpfr_ptr) 0);
+          }  /* prec */
+      }  /* i */
   mpfr_clears (x, y, (mpfr_ptr) 0);
   set_emin (emin);
 }
@@ -242,6 +247,76 @@
   large (MPFR_EMAX_MAX);
 }
 
+/* Cases where the function overflows on n = 0 when rounding is like
+   away from zero. */
+static void
+overflow0 (mpfr_exp_t emax)
+{
+  mpfr_exp_t old_emax;
+  mpfr_t x, y1, y2;
+  int neg, r, op;
+  static char *sop[4] = { "mul_2ui", "mul_2si", "div_2ui", "div_2si" };
+
+  old_emax = mpfr_get_emax ();
+  set_emax (emax);
+
+  mpfr_init2 (x, 8);
+  mpfr_inits2 (6, y1, y2, (mpfr_ptr) 0);
+
+  mpfr_set_inf (x, 1);
+  mpfr_nextbelow (x);
+
+  for (neg = 0; neg <= 1; neg++)
+    {
+      RND_LOOP (r)
+        {
+          int inex1, inex2;
+          unsigned int flags1, flags2;
+
+          /* Even if there isn't an overflow (rounding ~ toward zero),
+             the result is the same as the one of an overflow. */
+          inex1 = mpfr_overflow (y1, (mpfr_rnd_t) r, neg ? -1 : 1);
+          flags1 = MPFR_FLAGS_INEXACT;
+          if (mpfr_inf_p (y1))
+            flags1 |= MPFR_FLAGS_OVERFLOW;
+          for (op = 0; op < 4; op++)
+            {
+              mpfr_clear_flags ();
+              inex2 =
+                op == 0 ? mpfr_mul_2ui (y2, x, 0, (mpfr_rnd_t) r) :
+                op == 1 ? mpfr_mul_2si (y2, x, 0, (mpfr_rnd_t) r) :
+                op == 2 ? mpfr_div_2ui (y2, x, 0, (mpfr_rnd_t) r) :
+                op == 3 ? mpfr_div_2si (y2, x, 0, (mpfr_rnd_t) r) :
+                (MPFR_ASSERTN (0), 0);
+              flags2 = __gmpfr_flags;
+              if (!(mpfr_equal_p (y1, y2) &&
+                    SAME_SIGN (inex1, inex2) &&
+                    flags1 == flags2))
+                {
+                  printf ("Error in overflow0 for %s, mpfr_%s, emax = %"
+                          MPFR_EXP_FSPEC "d,\nx = ",
+                          mpfr_print_rnd_mode ((mpfr_rnd_t) r), sop[op],
+                          (mpfr_eexp_t) emax);
+                  mpfr_dump (x);
+                  printf ("Expected ");
+                  mpfr_dump (y1);
+                  printf ("  with inex = %d, flags =", inex1);
+                  flags_out (flags1);
+                  printf ("Got      ");
+                  mpfr_dump (y2);
+                  printf ("  with inex = %d, flags =", inex2);
+                  flags_out (flags2);
+                  exit (1);
+                }
+            }
+        }
+      mpfr_neg (x, x, MPFR_RNDN);
+    }
+
+  mpfr_clears (x, y1, y2, (mpfr_ptr) 0);
+  set_emax (old_emax);
+}
+
 int
 main (int argc, char *argv[])
 {
@@ -334,6 +409,11 @@
   underflow0 ();
   large0 ();
 
+  if (mpfr_get_emax () != MPFR_EMAX_MAX)
+    overflow0 (mpfr_get_emax ());
+  overflow0 (MPFR_EMAX_MAX);
+  overflow0 (-1);
+
   tests_end_mpfr ();
   return 0;
 }
--- a/tests/troot.c	2015-06-19 21:55:10.000000000 +0200
+++ b/tests/troot.c	2016-02-16 14:55:39.875392372 +0100
@@ -25,6 +25,19 @@
 
 #include "mpfr-test.h"
 
+#define DEFN(N)                                                         \
+  static int root##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)        \
+  { return mpfr_root (y, x, N, rnd); }                                  \
+  static int pow##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)         \
+  { return mpfr_pow_ui (y, x, N, rnd); }
+
+DEFN(2)
+DEFN(3)
+DEFN(4)
+DEFN(5)
+DEFN(17)
+DEFN(120)
+
 static void
 special (void)
 {
@@ -52,7 +65,7 @@
       exit (1);
     }
 
-  /* root(-Inf, 17) =  -Inf */
+  /* root(-Inf, 17) = -Inf */
   mpfr_set_inf (x, -1);
   mpfr_root (y, x, 17, MPFR_RNDN);
   if (!mpfr_inf_p (y) || mpfr_sgn (y) > 0)
@@ -69,7 +82,7 @@
       exit (1);
     }
 
-  /* root(+/-0) =  +/-0 */
+  /* root(+/-0, k) = +/-0 for k > 0 */
   mpfr_set_ui (x, 0, MPFR_RNDN);
   mpfr_root (y, x, 17, MPFR_RNDN);
   if (mpfr_cmp_ui (y, 0) || mpfr_sgn (y) < 0)
@@ -190,64 +203,39 @@
   i = mpfr_root (y, x, 1, MPFR_RNDN);
   if (mpfr_cmp_ui (x, 17) || i != 0)
     {
-      printf ("Error in root (17^(1/1))\n");
+      printf ("Error in root for 17^(1/1)\n");
       exit (1);
     }
 
-#if 0
-  /* Check for k == 0:
-     For 0 <= x < 1 => +0.
-     For x = 1      => 1.
-     For x > 1,     => +Inf.
-     For x < 0      => NaN.   */
-  i = mpfr_root (y, x, 0, MPFR_RNDN);
-  if (!MPFR_IS_INF (y) || !MPFR_IS_POS (y) || i != 0)
-    {
-      printf ("Error in root 17^(1/0)\n");
-      exit (1);
-    }
-  mpfr_set_ui (x, 1, MPFR_RNDN);
-  i = mpfr_root (y, x, 0, MPFR_RNDN);
-  if (mpfr_cmp_ui (y, 1) || i != 0)
-    {
-      printf ("Error in root 1^(1/0)\n");
-      exit (1);
-    }
   mpfr_set_ui (x, 0, MPFR_RNDN);
   i = mpfr_root (y, x, 0, MPFR_RNDN);
-  if (!MPFR_IS_ZERO (y) || !MPFR_IS_POS (y) || i != 0)
-    {
-      printf ("Error in root 0+^(1/0)\n");
-      exit (1);
-    }
-  MPFR_CHANGE_SIGN (x);
-  i = mpfr_root (y, x, 0, MPFR_RNDN);
-  if (!MPFR_IS_ZERO (y) || !MPFR_IS_POS (y) || i != 0)
+  if (!MPFR_IS_NAN (y) || i != 0)
     {
-      printf ("Error in root 0-^(1/0)\n");
+      printf ("Error in root for (+0)^(1/0)\n");
       exit (1);
     }
-  mpfr_set_ui_2exp (x, 17, -5, MPFR_RNDD);
+  mpfr_neg (x, x, MPFR_RNDN);
   i = mpfr_root (y, x, 0, MPFR_RNDN);
-  if (!MPFR_IS_ZERO (y) || !MPFR_IS_POS (y) || i != 0)
+  if (!MPFR_IS_NAN (y) || i != 0)
     {
-      printf ("Error in root (17/2^5)^(1/0)\n");
+      printf ("Error in root for (-0)^(1/0)\n");
       exit (1);
     }
-#endif
-  mpfr_set_ui (x, 0, MPFR_RNDN);
+
+  mpfr_set_ui (x, 1, MPFR_RNDN);
   i = mpfr_root (y, x, 0, MPFR_RNDN);
   if (!MPFR_IS_NAN (y) || i != 0)
     {
-      printf ("Error in root 0+^(1/0)\n");
+      printf ("Error in root for 1^(1/0)\n");
       exit (1);
     }
+
   /* Check for k==2 */
   mpfr_set_si (x, -17, MPFR_RNDD);
   i = mpfr_root (y, x, 2, MPFR_RNDN);
   if (!MPFR_IS_NAN (y) || i != 0)
     {
-      printf ("Error in root (-17)^(1/2)\n");
+      printf ("Error in root for (-17)^(1/2)\n");
       exit (1);
     }
 
@@ -255,11 +243,168 @@
   mpfr_clear (y);
 }
 
+/* https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=812779
+ * https://bugzilla.gnome.org/show_bug.cgi?id=756960
+ * is a GNOME Calculator bug (mpfr_root applied on a negative integer,
+ * which is converted to an unsigned integer), but the strange result
+ * is also due to a bug in MPFR.
+ */
+static void
+bigint (void)
+{
+  mpfr_t x, y;
+
+  mpfr_inits2 (64, x, y, (mpfr_ptr) 0);
+
+  mpfr_set_ui (x, 10, MPFR_RNDN);
+  if (sizeof (unsigned long) * CHAR_BIT == 64)
+    {
+      mpfr_root (x, x, ULONG_MAX, MPFR_RNDN);
+      mpfr_set_ui_2exp (y, 1, -63, MPFR_RNDN);
+      mpfr_add_ui (y, y, 1, MPFR_RNDN);
+      if (! mpfr_equal_p (x, y))
+        {
+          printf ("Error in bigint for ULONG_MAX\n");
+          printf ("Expected ");
+          mpfr_dump (y);
+          printf ("Got      ");
+          mpfr_dump (x);
+          exit (1);
+        }
+    }
+
+  mpfr_set_ui (x, 10, MPFR_RNDN);
+  mpfr_root (x, x, 1234567890, MPFR_RNDN);
+  mpfr_set_str_binary (y,
+    "1.00000000000000000000000000001000000000101011000101000110010001");
+  if (! mpfr_equal_p (x, y))
+    {
+      printf ("Error in bigint for 1234567890\n");
+      printf ("Expected ");
+      mpfr_dump (y);
+      printf ("Got      ");
+      mpfr_dump (x);
+      exit (1);
+    }
+
+  mpfr_clears (x, y, (mpfr_ptr) 0);
+}
+
 #define TEST_FUNCTION mpfr_root
 #define INTEGER_TYPE unsigned long
-#define INT_RAND_FUNCTION() (INTEGER_TYPE) (randlimb () % 3 +2)
+#define INT_RAND_FUNCTION() \
+  (INTEGER_TYPE) (randlimb () & 1 ? randlimb () : randlimb () % 3 + 2)
 #include "tgeneric_ui.c"
 
+static void
+exact_powers (unsigned long bmax, unsigned long kmax)
+{
+  long b, k;
+  mpz_t z;
+  mpfr_t x, y;
+  int inex, neg;
+
+  mpz_init (z);
+  for (b = 2; b <= bmax; b++)
+    for (k = 1; k <= kmax; k++)
+      {
+        mpz_ui_pow_ui (z, b, k);
+        mpfr_init2 (x, mpz_sizeinbase (z, 2));
+        mpfr_set_ui (x, b, MPFR_RNDN);
+        mpfr_pow_ui (x, x, k, MPFR_RNDN);
+        mpz_set_ui (z, b);
+        mpfr_init2 (y, mpz_sizeinbase (z, 2));
+        for (neg = 0; neg <= 1; neg++)
+          {
+            inex = mpfr_root (y, x, k, MPFR_RNDN);
+            if (inex != 0)
+              {
+                printf ("Error in exact_powers, b=%ld, k=%ld\n", b, k);
+                printf ("Expected inex=0, got %d\n", inex);
+                exit (1);
+              }
+            if (neg && (k & 1) == 0)
+              {
+                if (!MPFR_IS_NAN (y))
+                  {
+                    printf ("Error in exact_powers, b=%ld, k=%ld\n", b, k);
+                    printf ("Expected y=NaN\n");
+                    printf ("Got      ");
+                    mpfr_out_str (stdout, 10, 0, y, MPFR_RNDN);
+                    printf ("\n");
+                    exit (1);
+                  }
+              }
+            else if (MPFR_IS_NAN (y) || mpfr_cmp_si (y, b) != 0)
+              {
+                printf ("Error in exact_powers, b=%ld, k=%ld\n", b, k);
+                printf ("Expected y=%ld\n", b);
+                printf ("Got      ");
+                mpfr_out_str (stdout, 10, 0, y, MPFR_RNDN);
+                printf ("\n");
+                exit (1);
+              }
+            mpfr_neg (x, x, MPFR_RNDN);
+            b = -b;
+          }
+        mpfr_clear (x);
+        mpfr_clear (y);
+      }
+  mpz_clear (z);
+}
+
+/* Compare root(x,2^h) with pow(x,2^(-h)). */
+static void
+cmp_pow (void)
+{
+  mpfr_t x, y1, y2;
+  int h;
+
+  mpfr_inits2 (128, x, y1, y2, (mpfr_ptr) 0);
+
+  for (h = 1; h < sizeof (unsigned long) * CHAR_BIT; h++)
+    {
+      unsigned long k = (unsigned long) 1 << h;
+      int i;
+
+      for (i = 0; i < 10; i++)
+        {
+          mpfr_rnd_t rnd;
+          unsigned int flags1, flags2;
+          int inex1, inex2;
+
+          tests_default_random (x, 0, __gmpfr_emin, __gmpfr_emax, 1);
+          rnd = RND_RAND ();
+          mpfr_set_ui_2exp (y1, 1, -h, MPFR_RNDN);
+          mpfr_clear_flags ();
+          inex1 = mpfr_pow (y1, x, y1, rnd);
+          flags1 = __gmpfr_flags;
+          mpfr_clear_flags ();
+          inex2 = mpfr_root (y2, x, k, rnd);
+          flags2 = __gmpfr_flags;
+          if (!(mpfr_equal_p (y1, y2) && SAME_SIGN (inex1, inex2) &&
+                flags1 == flags2))
+            {
+              printf ("Error in cmp_pow on h=%d, i=%d, rnd=%s\n",
+                      h, i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
+              printf ("x = ");
+              mpfr_dump (x);
+              printf ("pow  = ");
+              mpfr_dump (y1);
+              printf ("with inex = %d, flags =", inex1);
+              flags_out (flags1);
+              printf ("root = ");
+              mpfr_dump (y2);
+              printf ("with inex = %d, flags =", inex2);
+              flags_out (flags2);
+              exit (1);
+            }
+        }
+    }
+
+  mpfr_clears (x, y1, y2, (mpfr_ptr) 0);
+}
+
 int
 main (void)
 {
@@ -270,7 +415,10 @@
 
   tests_start_mpfr ();
 
+  exact_powers (3, 1000);
   special ();
+  bigint ();
+  cmp_pow ();
 
   mpfr_init (x);
 
@@ -329,6 +477,13 @@
 
   test_generic_ui (2, 200, 30);
 
+  bad_cases (root2, pow2, "mpfr_root[2]", 8, -256, 255, 4, 128, 800, 40);
+  bad_cases (root3, pow3, "mpfr_root[3]", 8, -256, 255, 4, 128, 800, 40);
+  bad_cases (root4, pow4, "mpfr_root[4]", 8, -256, 255, 4, 128, 800, 40);
+  bad_cases (root5, pow5, "mpfr_root[5]", 8, -256, 255, 4, 128, 800, 40);
+  bad_cases (root17, pow17, "mpfr_root[17]", 8, -256, 255, 4, 128, 800, 40);
+  bad_cases (root120, pow120, "mpfr_root[120]", 8, -256, 255, 4, 128, 800, 40);
+
   tests_end_mpfr ();
   return 0;
 }
--- a/tests/tsqrt.c	2015-06-19 21:55:10.000000000 +0200
+++ b/tests/tsqrt.c	2016-02-16 14:55:39.839391671 +0100
@@ -569,6 +569,35 @@
   mpfr_clear (y);
 }
 
+/* Bug reported by Fredrik Johansson, occurring when:
+   - the precision of the result is a multiple of the number of bits
+     per word (GMP_NUMB_BITS),
+   - the rounding mode is to nearest (MPFR_RNDN),
+   - internally, the result has to be rounded up to a power of 2.
+*/
+static void
+bug20160120 (void)
+{
+  mpfr_t x, y;
+
+  mpfr_init2 (x, 4 * GMP_NUMB_BITS);
+  mpfr_init2 (y, GMP_NUMB_BITS);
+
+  mpfr_set_ui (x, 1, MPFR_RNDN);
+  mpfr_nextbelow (x);
+  mpfr_sqrt (y, x, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_check (y));
+  MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0);
+
+  mpfr_set_prec (y, 2 * GMP_NUMB_BITS);
+  mpfr_sqrt (y, x, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_check (y));
+  MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0);
+
+  mpfr_clear(x);
+  mpfr_clear(y);
+}
+
 #define TEST_FUNCTION test_sqrt
 #define TEST_RANDOM_POS 8
 #include "tgeneric.c"
@@ -704,6 +733,8 @@
   data_check ("data/sqrt", mpfr_sqrt, "mpfr_sqrt");
   bad_cases (mpfr_sqrt, mpfr_sqr, "mpfr_sqrt", 8, -256, 255, 4, 128, 800, 50);
 
+  bug20160120 ();
+
   tests_end_mpfr ();
   return 0;
 }
--- a/tests/tzeta.c	2015-06-19 21:55:10.000000000 +0200
+++ b/tests/tzeta.c	2016-02-16 14:55:39.834391573 +0100
@@ -394,6 +394,27 @@
   mpfr_nextabove (s);
   MPFR_ASSERTN (mpfr_equal_p (z, s) && inex > 0);
 
+  /* bug reported by Fredrik Johansson on 19 Jan 2016 */
+  mpfr_set_prec (s, 536);
+  mpfr_set_ui_2exp (s, 1, -424, MPFR_RNDN);
+  mpfr_sub_ui (s, s, 128, MPFR_RNDN);  /* -128 + 2^(-424) */
+  for (prec = 6; prec <= 536; prec += 8) /* should go through 318 */
+    {
+      mpfr_set_prec (z, prec);
+      mpfr_zeta (z, s, MPFR_RNDD);
+      mpfr_set_prec (y, prec + 10);
+      mpfr_zeta (y, s, MPFR_RNDD);
+      mpfr_prec_round (y, prec, MPFR_RNDD);
+      if (! mpfr_equal_p (z, y))
+        {
+          printf ("mpfr_zeta fails near -128 for inprec=%lu outprec=%lu\n",
+                  (unsigned long) mpfr_get_prec (s), (unsigned long) prec);
+          printf ("expected "); mpfr_dump (y);
+          printf ("got      "); mpfr_dump (z);
+          exit (1);
+        }
+    }
+
   mpfr_clear (s);
   mpfr_clear (y);
   mpfr_clear (z);
--- a/VERSION	2015-06-19 21:55:09.000000000 +0200
+++ b/VERSION	2016-02-16 14:55:39.870392274 +0100
@@ -1 +1 @@
-3.1.3
+3.1.3-p12
