added mpfr_atanu and mpfr_atanpi

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

7 files changed, 359 insertions, 9 deletions

diff --git a/NEWS b/NEWS
index 799f77b7d..d2863d3db 100644
--- a/NEWS
+++ b/NEWS
@@ -22,9 +22,10 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
 
 Changes from versions 4.1.* to version 4.2.0:
 - The "" release.
-- New functions mpfr_cosu, mpfr_sinu, mpfr_tanu, mpfr_acosu and mpfr_asinu.
-- New functions mpfr_cospi, mpfr_sinpi, mpfr_tanpi, mpfr_acospi and
-  mpfr_asinpi.
+- New functions mpfr_cosu, mpfr_sinu, mpfr_tanu, mpfr_acosu, mpfr_asinu and
+  mpfr_atanu.
+- New functions mpfr_cospi, mpfr_sinpi, mpfr_tanpi, mpfr_acospi, mpfr_asinpi
+  and mpfr_atanpi.
 - New function mpfr_fmod_ui.
 - Bug fixes.
   In particular, for the formatted output functions (mpfr_printf, etc.),
diff --git a/doc/mpfr.texi b/doc/mpfr.texi
index e2f88f3af..a669adfef 100644
--- a/doc/mpfr.texi
+++ b/doc/mpfr.texi
@@ -2329,17 +2329,19 @@ small precision of @var{rop}.
 
 @deftypefun int mpfr_acosu (mpfr_t @var{rop}, mpfr_t @var{op}, unsigned long int @var{u}, mpfr_rnd_t @var{rnd})
 @deftypefunx int mpfr_asinu (mpfr_t @var{rop}, mpfr_t @var{op}, unsigned long int @var{u}, mpfr_rnd_t @var{rnd})
+@deftypefunx int mpfr_atanu (mpfr_t @var{rop}, mpfr_t @var{op}, unsigned long int @var{u}, mpfr_rnd_t @var{rnd})
 Set @var{rop} to @var{a} multiplied by @var{u} and divided by
-@m{2\pi,2@tie{}Pi}, where @var{a} is the arc-cosine (resp.@: arc-sine)
-of @var{op}.
+@m{2\pi,2@tie{}Pi}, where @var{a} is the arc-cosine (resp.@: arc-sine and
+arc-tangent) of @var{op}.
 For example, if @var{u} equals 360, @code{mpfr_acosu} yields the arc-cosine in
 degrees.
 @end deftypefun
 
 @deftypefun int mpfr_acospi (mpfr_t @var{rop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
 @deftypefunx int mpfr_asinpi (mpfr_t @var{rop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
-Set @var{rop} to @code{acos(op)} (resp.@: @code{asin(op)}) divided
-by @m{\pi,Pi}.
+@deftypefunx int mpfr_atnnpi (mpfr_t @var{rop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
+Set @var{rop} to @code{acos(op)} (resp.@: @code{asin(op)} and @code{atnn(op)})
+divided by @m{\pi,Pi}.
 @end deftypefun
 
 @deftypefun int mpfr_atan2 (mpfr_t @var{rop}, mpfr_t @var{y}, mpfr_t @var{x}, mpfr_rnd_t @var{rnd})
@@ -4123,6 +4125,8 @@ that were added after MPFR@tie{}2.2, and in which MPFR version.
 
 @item @code{mpfr_asprintf} in MPFR@tie{}2.4.
 
+@item @code{mpfr_atanpi} and @code{mpfr_atanu} in MPFR@tie{}4.2.0.
+
 @item @code{mpfr_beta} in MPFR@tie{}4.0 (incomplete, experimental).
 
 @item @code{mpfr_buildopt_decimal_p} in MPFR@tie{}3.0.
diff --git a/src/Makefile.am b/src/Makefile.am
index 7f74f7bb6..0c7f214f7 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -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
+acosu.c asinu.c atanu.c
 
 nodist_libmpfr_la_SOURCES = $(BUILT_SOURCES)
 
diff --git a/src/atanu.c b/src/atanu.c
new file mode 100644
index 000000000..3c4089b25
--- /dev/null
+++ b/src/atanu.c
@@ -0,0 +1,187 @@
+/* mpfr_atanu  -- atanu(x)  = atan(x)*u/(2*pi)
+   mpfr_atanpi -- atanpi(x) = atan(x)/pi
+
+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. */
+
+#define MPFR_NEED_LONGLONG_H
+#include "mpfr-impl.h"
+
+/* put in y the correctly rounded value of atan(x)*u/(2*pi) */
+int
+mpfr_atanu (mpfr_ptr y, mpfr_srcptr x, unsigned long u, mpfr_rnd_t rnd_mode)
+{
+  mpfr_t tmp, pi;
+  mpfr_prec_t prec;
+  mpfr_exp_t expx;
+  int inex;
+  MPFR_SAVE_EXPO_DECL (expo);
+  MPFR_ZIV_DECL (loop);
+
+  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));
+
+  /* Singular cases */
+  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
+    {
+      if (MPFR_IS_NAN (x))
+        {
+          MPFR_SET_NAN (y);
+          MPFR_RET_NAN;
+        }
+      else if (MPFR_IS_INF (x))
+        {
+          /* atanu(+Inf,u) = u/4, atanu(-Inf,u) = -u/4 */
+          if (MPFR_IS_POS (x))
+            return mpfr_set_ui_2exp (y, u, -2, rnd_mode);
+          else
+            {
+              inex = mpfr_set_ui_2exp (y, u, -2, MPFR_INVERT_RND (rnd_mode));
+              MPFR_CHANGE_SIGN (y);
+              return -inex;
+            }
+        }
+      else /* necessarily x=0 */
+        {
+          MPFR_ASSERTD(MPFR_IS_ZERO(x));
+          /* atan(0)=0 with same sign, even when u=0 to ensure
+             atanu(-x,u) = -atanu(x,u) */
+          MPFR_SET_ZERO (y);
+          MPFR_SET_SAME_SIGN (y, x);
+          MPFR_RET (0); /* exact result */
+        }
+    }
+
+  if (u == 0) /* return 0 with sign of x, which is coherent with case x=0 */
+    {
+      MPFR_SET_ZERO (y);
+      MPFR_SET_POS (y);
+      MPFR_RET (0);
+    }
+
+  if (mpfr_cmpabs_ui (x, 1) == 0)
+    {
+      /* |x| = 1: atanu(1,u) = u/8, atanu(-1,u)=-u/8 */
+      /* we can't use mpfr_set_si_2exp with -u since -u might not be
+         representable as long */
+      if (MPFR_SIGN(x) > 0)
+        return mpfr_set_ui_2exp (y, u, -3, rnd_mode);
+      else
+        {
+          inex = mpfr_set_ui_2exp (y, u, -3, MPFR_INVERT_RND(rnd_mode));
+          MPFR_CHANGE_SIGN(y);
+          return -inex;
+        }
+    }
+
+  /* For x>=1, we have pi/2-1/x < atan(x) < pi/2, thus
+     u/4-u/(2*pi*x) < atanu(x,u) < u/4, and the relative difference between
+     atanu(x,u) and u/4 is less than 2/(pi*x) < 1/x <= 2^(1-EXP(x)).
+     If the relative difference is <= 2^(-prec-2), then the difference
+     between atanu(x,u) and u/4 is <= 1/4*ulp(u/4) <= 1/2*ulp(RN(u/4)).
+     We also require x >= 2^64, which implies x > 2*u/pi, so that
+     (u-1)/4 < u/4-u/(2*pi*x) < u/4. */
+  expx = MPFR_GET_EXP(x);
+  if (expx >= 65 && expx - 1 >= MPFR_PREC(y) + 2)
+    {
+      prec = (MPFR_PREC(y) <= 63) ? 65 : MPFR_PREC(y) + 2;
+      /* now prec > 64 and prec > MPFR_PREC(y)+1 */
+      mpfr_init2 (tmp, prec);
+      inex = mpfr_set_ui (tmp, u, MPFR_RNDN); /* exact since prec >= 64 */
+      MPFR_ASSERTD(inex == 0);
+      mpfr_nextbelow (tmp);
+      /* Since prec >= 65, the last significant bit of tmp is 1, and since
+         prec > PREC(y), tmp is not representable in the target precision,
+         which ensures we will get a correct ternary value below. */
+      MPFR_ASSERTD(mpfr_min_prec(tmp) > MPFR_PREC(y));
+      if (MPFR_SIGN(x) < 0)
+        MPFR_CHANGE_SIGN(tmp);
+      /* since prec >= PREC(y)+2, the rounding of tmp is correct */
+      inex = mpfr_div_2ui (y, tmp, 2, rnd_mode);
+      mpfr_clear (tmp);
+      return inex;
+    }
+
+  prec = MPFR_PREC (y);
+
+  MPFR_SAVE_EXPO_MARK (expo);
+
+  prec += MPFR_INT_CEIL_LOG2(prec) + 10;
+
+  mpfr_init2 (tmp, prec);
+  mpfr_init2 (pi, prec);
+
+  MPFR_ZIV_INIT (loop, prec);
+  for (;;)
+    {
+      /* In the error analysis below, each thetax denotes a variable such that
+         |thetax| <= 2^(1-prec) */
+      mpfr_atan (tmp, x, MPFR_RNDA);
+      /* tmp = atan(x) * (1 + theta1), and tmp cannot be zero since we rounded
+         away from zero, and the case x=0 was treated before */
+      /* first multiply by u to avoid underflow issues */
+      mpfr_mul_ui (tmp, tmp, u, MPFR_RNDA);
+      /* tmp = atan(x)*u * (1 + theta2)^2, and |tmp| >= 0.5*2^emin */
+      mpfr_const_pi (pi, MPFR_RNDZ); /* round toward zero since we we will
+                                        divide by pi, to round tmp away */
+      /* pi = Pi * (1 + theta3) */
+      mpfr_div (tmp, tmp, pi, MPFR_RNDA);
+      /* tmp = atan(x)*u/Pi * (1 + theta4)^4, with |tmp| > 0 */
+      /* since we rounded away from 0, if we get 0.5*2^emin here, it means
+         |atanu(x,u)| < 0.25*2^emin (pi is not exact) thus we have underflow */
+      if (MPFR_EXP(tmp) == __gmpfr_emin)
+        {
+          /* mpfr_underflow rounds away for RNDN */
+          mpfr_clear (tmp);
+          mpfr_clear (pi);
+          MPFR_SAVE_EXPO_FREE (expo);
+          return mpfr_underflow (y,
+                            (rnd_mode == MPFR_RNDN) ? MPFR_RNDZ : rnd_mode, 1);
+        }
+      mpfr_div_2ui (tmp, tmp, 1, MPFR_RNDA); /* exact */
+      /* tmp = atan(x)*u/(2*Pi) * (1 + theta4)^4 */
+      /* since |(1 + theta4)^4 - 1| <= 8*|theta4| for prec >= 3,
+         the relative error is less than 2^(4-prec) */
+      MPFR_ASSERTD(!MPFR_IS_ZERO(tmp));
+      if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - 4,
+                                       MPFR_PREC (y), rnd_mode)))
+        break;
+      MPFR_ZIV_NEXT (loop, prec);
+      mpfr_set_prec (tmp, prec);
+      mpfr_set_prec (pi, prec);
+    }
+  MPFR_ZIV_FREE (loop);
+
+  inex = mpfr_set (y, tmp, rnd_mode);
+  mpfr_clear (tmp);
+  mpfr_clear (pi);
+
+  MPFR_SAVE_EXPO_FREE (expo);
+  return mpfr_check_range (y, inex, rnd_mode);
+}
+
+int
+mpfr_atanpi (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
+{
+  return mpfr_atanu (y, x, 2, rnd_mode);
+}
diff --git a/src/mpfr.h b/src/mpfr.h
index 33151e9a4..c434856ec 100644
--- a/src/mpfr.h
+++ b/src/mpfr.h
@@ -728,8 +728,10 @@ __MPFR_DECLSPEC int mpfr_cosu (mpfr_ptr, mpfr_srcptr, unsigned long, mpfr_rnd_t)
 __MPFR_DECLSPEC int mpfr_tanu (mpfr_ptr, mpfr_srcptr, unsigned long, mpfr_rnd_t);
 __MPFR_DECLSPEC int mpfr_acosu (mpfr_ptr, mpfr_srcptr, unsigned long, mpfr_rnd_t);
 __MPFR_DECLSPEC int mpfr_asinu (mpfr_ptr, mpfr_srcptr, unsigned long, mpfr_rnd_t);
+__MPFR_DECLSPEC int mpfr_atanu (mpfr_ptr, mpfr_srcptr, unsigned long, mpfr_rnd_t);
 __MPFR_DECLSPEC int mpfr_acospi (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
 __MPFR_DECLSPEC int mpfr_asinpi (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
+__MPFR_DECLSPEC int mpfr_atanpi (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
 
 __MPFR_DECLSPEC int mpfr_sinpi (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
 __MPFR_DECLSPEC int mpfr_cospi (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
diff --git a/tests/Makefile.am b/tests/Makefile.am
index c0e4394c0..738abfa20 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -31,7 +31,8 @@ TESTS_NO_TVERSION = tabort_prec_max tassert tabort_defalloc1		\
      tabort_defalloc2 talloc tinternals tinits tisqrt tsgn tcheck	\
      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 taway tbeta tbuildopt \
+     talloc-cache tasin tasinh tasinu tatan tatanh tatanu taway tbeta 	\
+     tbuildopt 								\
      tcan_round tcbrt tcmp tcmp2 tcmp_d tcmp_ld tcmp_ui tcmpabs		\
      tcomparisons tconst_catalan tconst_euler tconst_log2 tconst_pi	\
      tcopysign tcos tcosh tcosu tcot tcoth tcsc tcsch td_div td_sub     \
diff --git a/tests/tatanu.c b/tests/tatanu.c
new file mode 100644
index 000000000..f69af4f5e
--- /dev/null
+++ b/tests/tatanu.c
@@ -0,0 +1,155 @@
+/* Test file for mpfr_atanu.
+
+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_atanu
+#define ULONG_ARG2
+#include "tgeneric.c"
+
+static void
+check_underflow (void)
+{
+  mpfr_t x, y;
+  mpfr_exp_t emin = mpfr_get_emin ();
+
+  set_emin (mpfr_get_emin_min ());
+
+  mpfr_init2 (x, MPFR_PREC_MIN);
+  mpfr_init2 (y, MPFR_PREC_MIN);
+  mpfr_set_ui_2exp (x, 1, mpfr_get_emin_min () - 1, MPFR_RNDN);
+  /* atanu(x,1) = atan(x)/(2*pi) will underflow */
+  mpfr_atanu (y, x, 1, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_zero_p (y) && mpfr_signbit (y) == 0);
+  mpfr_atanu (y, x, 1, MPFR_RNDZ);
+  MPFR_ASSERTN(mpfr_zero_p (y) && mpfr_signbit (y) == 0);
+  mpfr_atanu (y, x, 1, MPFR_RNDU);
+  MPFR_ASSERTN(mpfr_cmp_ui_2exp (y, 1, mpfr_get_emin_min () - 1) == 0);
+  mpfr_clear (x);
+  mpfr_clear (y);
+
+  set_emin (emin);
+}
+
+int
+main (void)
+{
+  mpfr_t x, y;
+  int r, inex;
+  unsigned long u;
+
+  tests_start_mpfr ();
+
+  check_underflow ();
+
+  mpfr_init (x);
+  mpfr_init (y);
+
+  /* check singular cases */
+  MPFR_SET_NAN(x);
+  mpfr_atanu (y, x, 1, MPFR_RNDN);
+  if (mpfr_nan_p (y) == 0)
+    {
+      printf ("Error: atanu (NaN, 1) != NaN\n");
+      exit (1);
+    }
+
+  mpfr_set_inf (x, 1);
+  mpfr_atanu (y, x, 1, MPFR_RNDN);
+  if (mpfr_cmp_ui_2exp (y, 1, -2) != 0)
+    {
+      printf ("Error: atanu (+Inf, 1) != 1/4\n");
+      printf ("got: "); mpfr_dump (y);
+      exit (1);
+    }
+
+  mpfr_set_inf (x, -1);
+  mpfr_atanu (y, x, 1, MPFR_RNDN);
+  if (mpfr_cmp_si_2exp (y, -1, -2) != 0)
+    {
+      printf ("Error: atanu (-Inf, 1) != -1/4\n");
+      exit (1);
+    }
+
+  /* atanu (+0,u) = +0 */
+  mpfr_set_ui (x, 0, MPFR_RNDN);
+  mpfr_atanu (y, x, 1, MPFR_RNDN);
+  if (MPFR_NOTZERO (y) || MPFR_IS_NEG (y))
+    {
+      printf ("Error: atanu(+0,1) != +0\n");
+      exit (1);
+    }
+
+  /* atanu (-0,u) = -0 */
+  mpfr_set_ui (x, 0, MPFR_RNDN);
+  mpfr_neg (x, x, MPFR_RNDN);
+  mpfr_atanu (y, x, 1, MPFR_RNDN);
+  if (MPFR_NOTZERO (y) || MPFR_IS_POS (y))
+    {
+      printf ("Error: atanu(-0,1) != -0\n");
+      exit (1);
+    }
+
+  /* atanu (1,u) = u/8 */
+  for (u = 1; u < 100; u++)
+     RND_LOOP (r)
+       {
+         mpfr_set_ui (x, 1, MPFR_RNDN); /* exact */
+         mpfr_atanu (y, x, u, (mpfr_rnd_t) r);
+         inex = mpfr_set_ui (x, u, MPFR_RNDN);
+         MPFR_ASSERTN(inex == 0);
+         inex = mpfr_div_ui (x, x, 8, (mpfr_rnd_t) r);
+         if ((r != MPFR_RNDF || inex == 0) && !mpfr_equal_p (x, y))
+           {
+             printf ("Error: atanu(1,u) != u/8 for u=%lu and rnd=%s\n",
+                     u, mpfr_print_rnd_mode ((mpfr_rnd_t) r));
+             printf ("got: "); mpfr_dump (y);
+             exit (1);
+           }
+       }
+
+  /* atanu (-1,u) = -u/8 */
+  for (u = 1; u < 100; u++)
+     RND_LOOP (r)
+       {
+         mpfr_set_si (x, -1, MPFR_RNDN); /* exact */
+         mpfr_atanu (y, x, u, (mpfr_rnd_t) r);
+         inex = mpfr_set_ui (x, u, MPFR_RNDN);
+         MPFR_ASSERTN(inex == 0);
+         inex = mpfr_div_si (x, x, -8, (mpfr_rnd_t) r);
+         if ((r != MPFR_RNDF || inex == 0) && !mpfr_equal_p (x, y))
+           {
+             printf ("Error: atanu(-1,u) != -u/8 for u=%lu and rnd=%s\n",
+                     u, mpfr_print_rnd_mode ((mpfr_rnd_t) r));
+             printf ("got: "); mpfr_dump (y);
+             exit (1);
+           }
+       }
+
+  test_generic (MPFR_PREC_MIN, 100, 100);
+
+  mpfr_clear (x);
+  mpfr_clear (y);
+
+  tests_end_mpfr ();
+  return 0;
+}