added log2p1 and compound (mpfr_compound is not finished yet)

git-svn-id: https://scm.gforge.inria.fr/anonscm/svn/mpfr/trunk@14331 280ebfd0-de03-0410-8827-d642c229c3f4

10 files changed, 534 insertions, 8 deletions

diff --git a/NEWS b/NEWS
index d2863d3db..4ef73579f 100644
--- a/NEWS
+++ b/NEWS
@@ -26,6 +26,7 @@ Changes from versions 4.1.* to version 4.2.0:
   mpfr_atanu.
 - New functions mpfr_cospi, mpfr_sinpi, mpfr_tanpi, mpfr_acospi, mpfr_asinpi
   and mpfr_atanpi.
+- New function log2p1.
 - New function mpfr_fmod_ui.
 - Bug fixes.
   In particular, for the formatted output functions (mpfr_printf, etc.),
diff --git a/TODO b/TODO
index 6328f586c..6693b2685 100644
--- a/TODO
+++ b/TODO
@@ -116,7 +116,7 @@ Table of contents:
   (page 392):
   atan2pi
   compoundn = (1+x)^n
-  exp10m1, exp2m1, log10p1, log2p1
+  exp10m1, exp2m1, log10p1
   pown (but there is mpfr_pow_si with n long instead of intmax_t)
   powr = exp(y*log(x)). The difference with pow = x^y is that when x < 0
   and y is an integer, powr returns NaN, also when x=1 and y=+/Inf,
diff --git a/doc/mpfr.texi b/doc/mpfr.texi
index 01a2db64d..67598d794 100644
--- a/doc/mpfr.texi
+++ b/doc/mpfr.texi
@@ -2210,8 +2210,9 @@ Set @var{rop} to @minus{}Inf if @var{op} is @pom{}0
 @end deftypefun
 
 @deftypefun int mpfr_log1p (mpfr_t @var{rop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
-Set @var{rop} to the logarithm of one plus @var{op},
-rounded in the direction @var{rnd}.
+@deftypefunx int mpfr_log2p1 (mpfr_t @var{rop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
+Set @var{rop} to the logarithm of one plus @var{op} (in radix two for
+@code{mpfr_log2p1}), rounded in the direction @var{rnd}.
 Set @var{rop} to @minus{}Inf if @var{op} is @minus{}1.
 @end deftypefun
 
@@ -4228,6 +4229,8 @@ use @code{mpfr_get_z_exp}.
 
 @item @code{mpfr_j0}, @code{mpfr_j1} and @code{mpfr_jn} in MPFR@tie{}2.3.
 
+@item @code{mpfr_log2p1} in MPFR@tie{}4.2.0.
+
 @item @code{mpfr_lgamma} in MPFR@tie{}2.3.
 
 @item @code{mpfr_li2} in MPFR@tie{}2.4.
diff --git a/src/Makefile.am b/src/Makefile.am
index 7242c8e1f..0f17dd80e 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -45,7 +45,7 @@ set_str.c set_str_raw.c set_ui.c set_z.c sqrt.c sqrt_ui.c sub.c sub1.c  \
 sub_ui.c rint.c ui_div.c ui_sub.c urandom.c urandomb.c get_z_2exp.c     \
 swap.c factorial.c cosh.c sinh.c tanh.c sinh_cosh.c acosh.c asinh.c     \
 atanh.c atan.c cmp2.c exp_2.c asin.c const_euler.c cos.c sin.c tan.c    \
-fma.c fms.c hypot.c log1p.c expm1.c log2.c log10.c ui_pow.c             \
+fma.c fms.c hypot.c log1p.c expm1.c log2.c log2p1.c log10.c ui_pow.c    \
 ui_pow_ui.c minmax.c dim.c signbit.c copysign.c setsign.c gmp_op.c      \
 init2.c acos.c sin_cos.c set_nan.c set_inf.c set_zero.c powerof2.c      \
 gamma.c set_ld.c get_ld.c cbrt.c volatile.c fits_s.h fits_sshort.c      \
@@ -68,7 +68,7 @@ random_deviate.h random_deviate.c erandom.c mpfr-mini-gmp.c             \
 mpfr-mini-gmp.h fmma.c log_ui.c gamma_inc.c ubf.c invert_limb.h         \
 invsqrt_limb.h beta.c odd_p.c get_q.c pool.c total_order.c set_d128.c   \
 get_d128.c nbits_ulong.c cmpabs_ui.c sinu.c cosu.c tanu.c fmod_ui.c     \
-acosu.c asinu.c atanu.c
+acosu.c asinu.c atanu.c compound.c
 
 nodist_libmpfr_la_SOURCES = $(BUILT_SOURCES)
 
diff --git a/src/compound.c b/src/compound.c
new file mode 100644
index 000000000..62c2acd92
--- /dev/null
+++ b/src/compound.c
@@ -0,0 +1,91 @@
+/* mpfr_compound --- compound(x,n) = (1+x)^n
+
+Copyright 2021 Free Software Foundation, Inc.
+Contributed by the AriC and Caramba projects, INRIA.
+
+This file is part of the GNU MPFR Library.
+
+The GNU MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+The GNU MPFR Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
+https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
+51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
+
+#include "mpfr-impl.h"
+
+/* put in z the correctly rounded value of (1+x)^n */
+int
+mpfr_compound (mpfr_ptr y, mpfr_srcptr x, long n, mpfr_rnd_t rnd_mode)
+{
+  int inex, compared;
+
+  MPFR_LOG_FUNC
+    (("x[%Pu]=%.*Rg u=%lu rnd=%d",
+      mpfr_get_prec(x), mpfr_log_prec, x,
+      u, rnd_mode),
+     ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, inex));
+
+  /* Special cases */
+  if (MPFR_IS_SINGULAR (x))
+    {
+      if (MPFR_IS_NAN (x) || (MPFR_IS_INF (x) && MPFR_SIGN (x) < 0))
+        {
+          MPFR_SET_NAN (y);
+          MPFR_RET_NAN;
+        }
+      else if (n == 0 || MPFR_IS_ZERO (x))
+        /* (1+Inf)^0 = 1 and (1+x)^0 = 1 */
+        return mpfr_set_ui (y, 1, rnd_mode);
+      else if (MPFR_IS_INF (x)) /* x = +Inf */
+        {
+          if (n < 0) /* (1+Inf)^n = +0 for n < 0 */
+            MPFR_SET_ZERO (y);
+          else /* n > 0: (1+Inf)^n = +Inf */
+            MPFR_SET_INF (y);
+          MPFR_SET_POS (y);
+          MPFR_RET (0); /* exact 0 or infinity */
+        }
+    }
+
+  /* (1+x)^n = NaN for x < -1 */
+  compared = mpfr_cmp_si (x, -1);
+  if (compared < 0)
+    {
+      MPFR_SET_NAN (y);
+      MPFR_RET_NAN;
+    }
+
+  /* compound(x,0) gives 1 for x >= 1 */
+  if (n == 0)
+    return mpfr_set_ui (y, 1, rnd_mode);
+
+  if (compared == 0)
+    {
+      if (n < 0)
+        {
+          /* compound(-1,n) = +Inf with divide-by-zero exception */
+          MPFR_SET_INF (y);
+          MPFR_SET_POS (y);
+          MPFR_SET_DIVBY0 ();
+          MPFR_RET (0);
+        }
+      else
+        {
+          /* compound(-1,n) = +0 */
+          MPFR_SET_ZERO (y);
+          MPFR_SET_POS (y);
+          MPFR_RET (0);
+        }
+    }
+
+  return mpfr_set_si (y, -17, rnd_mode);
+}
diff --git a/src/log2p1.c b/src/log2p1.c
new file mode 100644
index 000000000..a7d7a28e3
--- /dev/null
+++ b/src/log2p1.c
@@ -0,0 +1,157 @@
+/* mpfr_log2p1 -- Compute log2(1+x)
+
+Copyright 2001-2021 Free Software Foundation, Inc.
+Contributed by the AriC and Caramba projects, INRIA.
+
+This file is part of the GNU MPFR Library.
+
+The GNU MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+The GNU MPFR Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
+https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
+51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
+
+#define MPFR_NEED_LONGLONG_H /* needed for MPFR_INT_CEIL_LOG2 */
+#include "mpfr-impl.h"
+
+/* return non-zero if log2(1+x) is exactly representable in infinite precision,
+   and in such case the returned value is k such that 1+x = 2^k (the case k=0
+   cannot happen since we assume x<>0) */
+static mpfr_exp_t
+mpfr_log2p1_isexact (mpfr_srcptr x)
+{
+  /* log2(1+x) is exactly representable when 1+x is a power of two,
+     we thus simply compute 1+x with 1-bit precision and check whether
+     the addition is exact. This routine is called with extended exponent
+     range, thus no need to extend it. */
+  mpfr_t t;
+  int inex;
+  mpfr_init2 (t, 1);
+  inex = mpfr_add_ui (t, x, 1, MPFR_RNDZ);
+  mpfr_clear (t);
+  return (inex == 0) ? MPFR_GET_EXP(t) - 1 : 0;
+}
+
+/* in case x=2^k and we can decide of the correct rounding,
+   put the correctly-rounded value in y and return the corresponding
+   ternary value (which is necessarily non-zero),
+   otherwise return 0 */
+static int
+mpfr_log2p1_special (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
+{
+  mpfr_exp_t expx = MPFR_GET_EXP(x);
+  mpfr_exp_t k = expx - 1, expk;
+  mpfr_prec_t prec;
+  mpfr_t t;
+  int inex;
+
+  if (k <= 0 || mpfr_cmp_si_2exp (x, 1, k) != 0)
+    return 0;
+  /* k < log2(1+x) < k + 1/x/log(2) < k + 2/x */
+  expk = MPFR_INT_CEIL_LOG2(k); /* exponent of k */
+  /* 2/x < 2^(2-EXP(x)) thus if 2-EXP(x) < expk - PREC(y) - 1,
+     we have 2/x < 1/4*ulp(k) and we can decide the correct rounding */
+  if (2 - expx >= expk - MPFR_PREC(y) - 1)
+    return 0;
+  prec = (MPFR_PREC(y) <= 64) ? 66 : MPFR_PREC(y) + 2;
+  mpfr_init2 (t, prec);
+  mpfr_set_ui (t, k, MPFR_RNDZ); /* exact since prec >= 64 */
+  mpfr_nextabove (t);
+  /* now k < t < k + 2/x and round(t) = round(log2(1+x)) */
+  inex = mpfr_set (y, t, rnd_mode);
+  mpfr_clear (t);
+  /* Warning: for RNDF, the mpfr_set calls above might return 0 */
+  return (rnd_mode == MPFR_RNDF) ? 1 : inex;
+}
+
+/* The computation of log2p1 is done by log2p1(x) = log1p(x)/log(2) */
+int
+mpfr_log2p1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
+{
+  int comp, inexact, nloop;
+  mpfr_t t, lg2;
+  mpfr_prec_t Ny = MPFR_PREC(y), prec;
+  MPFR_ZIV_DECL (loop);
+  MPFR_SAVE_EXPO_DECL (expo);
+
+  MPFR_LOG_FUNC
+    (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
+     ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y,
+      inexact));
+
+  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
+    return mpfr_log1p (y, x, rnd_mode); /* same result for singular cases */
+
+  comp = mpfr_cmp_si (x, -1);
+  /* log2p1(x) is undefined for x < -1 */
+  if (MPFR_UNLIKELY(comp <= 0))
+    {
+      if (comp == 0)
+        /* x=0: log2p1(-1)=-inf (divide-by-zero exception) */
+        {
+          MPFR_SET_INF (y);
+          MPFR_SET_NEG (y);
+          MPFR_SET_DIVBY0 ();
+          MPFR_RET (0);
+        }
+      MPFR_SET_NAN (y);
+      MPFR_RET_NAN;
+    }
+
+  MPFR_SAVE_EXPO_MARK (expo);
+
+  prec = Ny + MPFR_INT_CEIL_LOG2 (Ny) + 6;
+
+  mpfr_init2 (t, prec);
+  mpfr_init2 (lg2, prec);
+
+  MPFR_ZIV_INIT (loop, prec);
+  for (nloop = 0; ; nloop++)
+    {
+      mpfr_log1p (t, x, MPFR_RNDN);
+      mpfr_const_log2 (lg2, MPFR_RNDN);
+      mpfr_div (t, t, lg2, MPFR_RNDN);
+      /* t = log2(1+x) * (1 + theta)^3 where |theta| < 2^-prec,
+         for prec >= 2 we have |(1 + theta)^3 - 1| < 4*theta. */
+      if (MPFR_LIKELY (MPFR_CAN_ROUND (t, prec - 2, Ny, rnd_mode)))
+        break;
+
+      if (nloop == 0)
+        {
+          /* check for exact cases */
+          mpfr_exp_t k = mpfr_log2p1_isexact (x);
+          if (k != 0) /* 1+x = 2^k */
+            {
+              inexact = mpfr_set_si (y, k, rnd_mode);
+              goto end;
+            }
+
+          /* if x = 2^k with huge k, Ziv's loop will fail */
+          inexact = mpfr_log2p1_special (y, x, rnd_mode);
+          if (inexact != 0)
+            goto end;
+        }
+
+      MPFR_ZIV_NEXT (loop, prec);
+      mpfr_set_prec (t, prec);
+      mpfr_set_prec (lg2, prec);
+    }
+  inexact = mpfr_set (y, t, rnd_mode);
+
+ end:
+  MPFR_ZIV_FREE (loop);
+  mpfr_clear (t);
+  mpfr_clear (lg2);
+
+  MPFR_SAVE_EXPO_FREE (expo);
+  return mpfr_check_range (y, inexact, rnd_mode);
+}
diff --git a/src/mpfr.h b/src/mpfr.h
index 10af1f4db..f7e39447a 100644
--- a/src/mpfr.h
+++ b/src/mpfr.h
@@ -555,6 +555,7 @@ __MPFR_DECLSPEC int mpfr_snprintf (char*, size_t, const char*, ...);
 
 __MPFR_DECLSPEC int mpfr_pow (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t);
 __MPFR_DECLSPEC int mpfr_pow_si (mpfr_ptr, mpfr_srcptr, long, mpfr_rnd_t);
+__MPFR_DECLSPEC int mpfr_compound (mpfr_ptr, mpfr_srcptr, long, mpfr_rnd_t);
 __MPFR_DECLSPEC int mpfr_pow_ui (mpfr_ptr, mpfr_srcptr, unsigned long,
                                  mpfr_rnd_t);
 __MPFR_DECLSPEC int mpfr_ui_pow_ui (mpfr_ptr, unsigned long, unsigned long,
@@ -612,6 +613,7 @@ __MPFR_DECLSPEC int mpfr_log (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
 __MPFR_DECLSPEC int mpfr_log2 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
 __MPFR_DECLSPEC int mpfr_log10 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
 __MPFR_DECLSPEC int mpfr_log1p (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
+__MPFR_DECLSPEC int mpfr_log2p1 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
 __MPFR_DECLSPEC int mpfr_log_ui (mpfr_ptr, unsigned long, mpfr_rnd_t);
 
 __MPFR_DECLSPEC int mpfr_exp (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 83524896e..57e4d4643 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -32,9 +32,9 @@ TESTS_NO_TVERSION = tabort_prec_max tassert tabort_defalloc1            \
      tisnan texceptions tset_exp tset mpf_compat mpfr_compat reuse      \
      tabs tacos tacosh tacosu tadd tadd1sp tadd_d tadd_ui tagm tai      \
      talloc-cache tasin tasinh tasinu tatan tatanh tatanu tatan2u taway \
-     tbeta tbuildopt                                                    \
-     tcan_round tcbrt tcmp tcmp2 tcmp_d tcmp_ld tcmp_ui tcmpabs         \
-     tcomparisons tconst_catalan tconst_euler tconst_log2 tconst_pi     \
+     tbeta tbuildopt tcan_round tcbrt tcmp tcmp2 tcmp_d tcmp_ld tcmp_ui \
+     tcmpabs tcomparisons tcompound tconst_catalan tconst_euler         \
+     tconst_log2 tconst_pi                                              \
      tcopysign tcos tcosh tcosu tcot tcoth tcsc tcsch td_div td_sub     \
      tdigamma tdim tdiv tdiv_d tdiv_ui tdot teint teq terandom          \
      terandom_chisq terf texp texp10 texp2 texpm1 tfactorial tfits tfma \
@@ -43,6 +43,7 @@ TESTS_NO_TVERSION = tabort_prec_max tassert tabort_defalloc1            \
      tget_set_d64 tget_set_d128 tget_sj tget_str tget_z tgmpop tgrandom \
      thyperbolic thypot tinp_str                                        \
      tj0 tj1 tjn tl2b tlgamma tli2 tlngamma tlog tlog10 tlog1p tlog2    \
+     tlog2p1                                                            \
      tlog_ui tmin_prec tminmax tmodf tmul tmul_2exp tmul_d tmul_ui      \
      tnext tnrandom tnrandom_chisq tout_str toutimpl tpow tpow3         \
      tpow_all tpow_z tprec_round tprintf trandom trandom_deviate        \
diff --git a/tests/tcompound.c b/tests/tcompound.c
new file mode 100644
index 000000000..4e661d519
--- /dev/null
+++ b/tests/tcompound.c
@@ -0,0 +1,173 @@
+/* Test file for mpfr_compound.
+
+Copyright 2021 Free Software Foundation, Inc.
+Contributed by the AriC and Caramba projects, INRIA.
+
+This file is part of the GNU MPFR Library.
+
+The GNU MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+The GNU MPFR Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
+https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
+51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
+
+#include "mpfr-test.h"
+
+#define TEST_FUNCTION mpfr_compound
+#define INTEGER_TYPE long
+#define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 1, RANDS)
+#define test_generic_ui test_generic_si
+#include "tgeneric_ui.c"
+
+/* Special cases from IEEE 754-2019 */
+static void
+check_ieee754 (void)
+{
+  mpfr_t x, y;
+  long i;
+  mpfr_prec_t prec = 2; /* we need at least 2 so that 3/4 is exact */
+
+  mpfr_init2 (x, prec);
+  mpfr_init2 (y, prec);
+
+  /* compound(x,n) = NaN for x < -1, and set invalid exception */
+  mpfr_clear_nanflag ();
+  mpfr_set_si (x, -2, MPFR_RNDN);
+  mpfr_compound (y, x, 17, MPFR_RNDN);
+  if (!mpfr_nan_p (y))
+    {
+      printf ("Error, compound(-2,17) should give NaN\n");
+      printf ("got "); mpfr_dump (y);
+      exit (1);
+    }
+  if (!mpfr_nanflag_p ())
+    {
+      printf ("Error, compound(-2,17) should raise invalid flag\n");
+      exit (1);
+    }
+
+  /* compound(x,0) = 1 for x >= -1 (or NaN), we choose 1 */
+  for (i = -1; i <= 2; i++)
+    {
+      if (i != 2)
+        mpfr_set_si (x, i, MPFR_RNDN);
+      else
+        mpfr_set_inf (x, 1);
+      mpfr_compound (y, x, 0, MPFR_RNDN);
+      if (mpfr_cmp_ui (y, 1) != 0)
+        {
+          mpfr_printf ("Error, compound(%Re,0) should give 1\n", x);
+          printf ("got "); mpfr_dump (y);
+          exit (1);
+        }
+    }
+
+  /* compound(-1,n) = +Inf for n < 0, and raise divide-by-zero flag */
+  mpfr_clear_divby0 ();
+  mpfr_set_si (x, -1, MPFR_RNDN);
+  mpfr_compound (y, x, -1, MPFR_RNDN);
+  if (!mpfr_inf_p (y) || MPFR_SIGN(y) < 0)
+    {
+      printf ("Error, compound(-1,-1) should give +Inf\n");
+      printf ("got "); mpfr_dump (y);
+      exit (1);
+    }
+  if (!mpfr_divby0_p ())
+    {
+      printf ("Error, compound(-1,-1) should raise divide-by-zero flag\n");
+      exit (1);
+    }
+
+  /* compound(-1,n) = +0 for n > 0 */
+  mpfr_set_si (x, -1, MPFR_RNDN);
+  mpfr_compound (y, x, 1, MPFR_RNDN);
+  if (!mpfr_zero_p (y) || MPFR_SIGN(y) < 0)
+    {
+      printf ("Error, compound(-1,1) should give +0\n");
+      printf ("got "); mpfr_dump (y);
+      exit (1);
+    }
+
+  /* compound(+/-0,n) = 1 */
+  for (i = -1; i <= 1; i++)
+    {
+      mpfr_set_zero (x, -1);
+      mpfr_compound (y, x, i, MPFR_RNDN);
+      if (mpfr_cmp_ui (y, 1) != 0)
+        {
+          mpfr_printf ("Error1, compound(%Re,%ld) should give 1\n", x, i);
+          printf ("got "); mpfr_dump (y);
+          exit (1);
+        }
+      mpfr_set_zero (x, +1);
+      mpfr_compound (y, x, i, MPFR_RNDN);
+      if (mpfr_cmp_ui (y, 1) != 0)
+        {
+          mpfr_printf ("Error, compound(%Re,%ld) should give 1\n", x, i);
+          printf ("got "); mpfr_dump (y);
+          exit (1);
+        }
+    }
+
+  /* compound(+Inf,n) = +Inf for n > 0 */
+  mpfr_set_inf (x, 1);
+  mpfr_compound (y, x, 1, MPFR_RNDN);
+  if (!mpfr_inf_p (y) || MPFR_SIGN(y) < 0)
+    {
+      printf ("Error, compound(+Inf,1) should give +Inf\n");
+      printf ("got "); mpfr_dump (y);
+      exit (1);
+    }
+
+  /* compound(+Inf,n) = +0 for n < 0 */
+  mpfr_set_inf (x, 1);
+  mpfr_compound (y, x, -1, MPFR_RNDN);
+  if (!mpfr_zero_p (y) || MPFR_SIGN(y) < 0)
+    {
+      printf ("Error, compound(+Inf,-1) should give +0\n");
+      printf ("got "); mpfr_dump (y);
+      exit (1);
+    }
+
+  /* compound(NaN,n) = NaN for n <> 0 */
+  mpfr_set_nan (x);
+  mpfr_compound (y, x, -1, MPFR_RNDN);
+  if (!mpfr_nan_p (y))
+    {
+      printf ("Error, compound(NaN,-1) should give NaN\n");
+      printf ("got "); mpfr_dump (y);
+      exit (1);
+    }
+  mpfr_compound (y, x, +1, MPFR_RNDN);
+  if (!mpfr_nan_p (y))
+    {
+      printf ("Error, compound(NaN,+1) should give NaN\n");
+      printf ("got "); mpfr_dump (y);
+      exit (1);
+    }
+
+  mpfr_clear (x);
+  mpfr_clear (y);
+}
+
+int
+main (void)
+{
+  tests_start_mpfr ();
+
+  check_ieee754 ();
+
+  test_generic_si (MPFR_PREC_MIN, 100, 100);
+
+  tests_end_mpfr ();
+  return 0;
+}
diff --git a/tests/tlog2p1.c b/tests/tlog2p1.c
new file mode 100644
index 000000000..1d703b37f
--- /dev/null
+++ b/tests/tlog2p1.c
@@ -0,0 +1,98 @@
+/* Test file for mpfr_log2p1.
+
+Copyright 2001-2021 Free Software Foundation, Inc.
+Contributed by the AriC and Caramba projects, INRIA.
+
+This file is part of the GNU MPFR Library.
+
+The GNU MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+The GNU MPFR Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
+https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
+51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
+
+#include "mpfr-test.h"
+
+#define TEST_FUNCTION mpfr_log2p1
+#define TEST_RANDOM_EMAX 80
+#include "tgeneric.c"
+
+static void
+special (void)
+{
+  mpfr_t x;
+  int inex;
+
+  mpfr_init2 (x, 32);
+
+  mpfr_set_nan (x);
+  mpfr_clear_flags ();
+  inex = mpfr_log2p1 (x, x, MPFR_RNDN);
+  MPFR_ASSERTN (mpfr_nan_p (x) && inex == 0);
+  MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_NAN);
+
+  mpfr_set_inf (x, -1);
+  mpfr_clear_flags ();
+  inex = mpfr_log2p1 (x, x, MPFR_RNDN);
+  MPFR_ASSERTN (mpfr_nan_p (x) && inex == 0);
+  MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_NAN);
+
+  mpfr_set_inf (x, 1);
+  mpfr_clear_flags ();
+  inex = mpfr_log2p1 (x, x, MPFR_RNDN);
+  MPFR_ASSERTN (mpfr_inf_p (x) && mpfr_sgn (x) > 0 && inex == 0);
+  MPFR_ASSERTN (__gmpfr_flags == 0);
+
+  mpfr_set_ui (x, 0, MPFR_RNDN);
+  mpfr_clear_flags ();
+  inex = mpfr_log2p1 (x, x, MPFR_RNDN);
+  MPFR_ASSERTN (mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS (x) && inex == 0);
+  MPFR_ASSERTN (__gmpfr_flags == 0);
+  mpfr_neg (x, x, MPFR_RNDN);
+  mpfr_clear_flags ();
+  inex = mpfr_log2p1 (x, x, MPFR_RNDN);
+  MPFR_ASSERTN (mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_NEG (x) && inex == 0);
+  MPFR_ASSERTN (__gmpfr_flags == 0);
+
+  mpfr_set_si (x, -1, MPFR_RNDN);
+  mpfr_clear_flags ();
+  inex = mpfr_log2p1 (x, x, MPFR_RNDN);
+  MPFR_ASSERTN (mpfr_inf_p (x) && mpfr_sgn (x) < 0 && inex == 0);
+  MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_DIVBY0);
+
+  mpfr_set_si (x, -2, MPFR_RNDN);
+  mpfr_clear_flags ();
+  inex = mpfr_log2p1 (x, x, MPFR_RNDN);
+  MPFR_ASSERTN (mpfr_nan_p (x) && inex == 0);
+  MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_NAN);
+
+  /* include one hard-coded test */
+  mpfr_set_ui (x, 17, MPFR_RNDN);
+  inex = mpfr_log2p1 (x, x, MPFR_RNDN);
+  MPFR_ASSERTN (mpfr_cmp_ui_2exp (x, 1119355719, -28) == 0);
+  MPFR_ASSERTN (inex < 0);
+
+  mpfr_clear (x);
+}
+
+int
+main (int argc, char *argv[])
+{
+  tests_start_mpfr ();
+
+  special ();
+
+  test_generic (MPFR_PREC_MIN, 100, 50);
+
+  tests_end_mpfr ();
+  return 0;
+}