added new function mpfr_cosu

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

7 files changed, 340 insertions, 14 deletions

diff --git a/NEWS b/NEWS
index 113f858a9..ec80ce0b0 100644
--- a/NEWS
+++ b/NEWS
@@ -22,7 +22,7 @@ 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 function mpfr_sinu.
+- New function mpfr_sinu and mpfr_cosu.
 - Bug fixes.
   In particular, for the formatted output functions (mpfr_printf, etc.),
   the case where the precision consists only of a period has been fixed
diff --git a/doc/mpfr.texi b/doc/mpfr.texi
index 0402ac537..afe2a25f9 100644
--- a/doc/mpfr.texi
+++ b/doc/mpfr.texi
@@ -2267,15 +2267,23 @@ Set @var{rop} to the cosine of @var{op}, sine of @var{op},
 tangent of @var{op}, rounded in the direction @var{rnd}.
 @end deftypefun
 
-@deftypefun int mpfr_sinu (mpfr_t @var{rop}, mpfr_t @var{op}, unsigned long @var{u}, mpfr_rnd_t @var{rnd})
-Set @var{rop} to the sine of @m{@var{op} \times 2\pi/u,@var{op} multiplied
-by 2@tie{}Pi and divided by @var{u}}. Thus for @var{u} equal to 360, one
-gets the sine for @var{op} in degrees. When @m{@var{op} \times 2/u,@var{op}
+@deftypefun int mpfr_cosu (mpfr_t @var{rop}, mpfr_t @var{op}, unsigned long @var{u}, mpfr_rnd_t @var{rnd})
+@deftypefunx int mpfr_sinu (mpfr_t @var{rop}, mpfr_t @var{op}, unsigned long @var{u}, mpfr_rnd_t @var{rnd})
+Set @var{rop} to the cosine (resp. sine) of @m{@var{op} \times 2\pi/u,@var{op} multiplied
+by 2@tie{}Pi and divided by @var{u}}. Thus if @var{u} equals 360, one
+gets the cosine (resp. sine) for @var{op} in degrees.
+For @code{mpfr_cosu},
+when @m{@var{op} \times 2/u,@var{op}
+multiplied by 2 and divided by @var{u}} is a half-integer, the result is +0,
+following IEEE@tie{}754-2019 (cosPi),
+so that the function is even.
+For @code{mpfr_sinu},
+when @m{@var{op} \times 2/u,@var{op}
 multiplied by 2 and divided by @var{u}} is an integer, the result is zero
 with the same sign as @var{op}, following IEEE@tie{}754-2019 (sinPi),
 so that the function is odd.
 
-Note: this function is experimental and its interface might change in future
+Note: these functions are experimental and their interface might change in future
 versions.
 @end deftypefun
 
@@ -4114,6 +4122,8 @@ but not documented,
 and with a slight difference in the semantics (when
 the second input operand is a NaN)@.
 
+@item @code{mpfr_cosu} in MPFR@tie{}4.2.0 (experimental).
+
 @item @code{mpfr_custom_get_significand} in MPFR@tie{}3.0.
 This function was named @code{mpfr_custom_get_mantissa} in previous
 versions; @code{mpfr_custom_get_mantissa} is still available via a
diff --git a/src/Makefile.am b/src/Makefile.am
index a96594693..10686c5f7 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -67,7 +67,7 @@ grandom.c fpif.c set_float128.c get_float128.c rndna.c nrandom.c        \
 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
+get_d128.c nbits_ulong.c cmpabs_ui.c sinu.c cosu.c
 
 nodist_libmpfr_la_SOURCES = $(BUILT_SOURCES)
 
diff --git a/src/cosu.c b/src/cosu.c
new file mode 100644
index 000000000..783553df5
--- /dev/null
+++ b/src/cosu.c
@@ -0,0 +1,141 @@
+/* mpfr_cosu -- cosu(x) = cos(2*pi*x/u)
+
+Copyright 2020 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"
+
+/* FIXME[VL]: Implement the range reduction in this function.
+   That's the whole point of cosu compared to cos. */
+
+/* put in y the corrected-rounded value of cos(2*pi*x/u) */
+int
+mpfr_cosu (mpfr_ptr y, mpfr_srcptr x, unsigned long u, mpfr_rnd_t rnd_mode)
+{
+  mpfr_prec_t precy, prec;
+  mpfr_exp_t expx, expt, err;
+  mpfr_t t;
+  int inexact = 0, nloops = 0, underflow = 0;
+  MPFR_ZIV_DECL (loop);
+  MPFR_SAVE_EXPO_DECL (expo);
+
+  if (u == 0 || MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
+    {
+      /* for u=0, return NaN */
+      if (u == 0 || MPFR_IS_NAN (x) || MPFR_IS_INF (x))
+        {
+          MPFR_SET_NAN (y);
+          MPFR_RET_NAN;
+        }
+      else /* x is zero: cos(0) = 1 */
+        {
+          MPFR_ASSERTD (MPFR_IS_ZERO (x));
+          return mpfr_set_ui (y, 1, rnd_mode);
+        }
+    }
+
+  MPFR_SAVE_EXPO_MARK (expo);
+
+  precy = MPFR_PREC (y);
+  expx = MPFR_GET_EXP (x);
+  /* for x large, since argument reduction is expensive, we want to avoid
+     any failure in Ziv's strategy, thus we take into account expx too */
+  prec = precy + MPFR_INT_CEIL_LOG2 (MAX(precy,expx)) + 8;
+  MPFR_ASSERTD(prec >= 2);
+  mpfr_init2 (t, prec);
+  MPFR_ZIV_INIT (loop, prec);
+  for (;;)
+    {
+      nloops ++;
+      /* We first compute an approximation t of 2*pi*x/u, then call cos(t).
+         If t = 2*pi*x/u + s, then |cos(t) - cos(2*pi*x/u)| <= |s|. */
+      mpfr_set_prec (t, prec);
+      mpfr_const_pi (t, MPFR_RNDN); /* t = pi * (1 + theta1) where
+                                       |theta1| <= 2^-prec */
+      mpfr_mul_2ui (t, t, 1, MPFR_RNDN); /* t = 2*pi * (1 + theta1) */
+      mpfr_mul (t, t, x, MPFR_RNDN);     /* t = 2*pi*x * (1 + theta2)^2 where
+                                            |theta2| <= 2^-prec */
+      mpfr_div_ui (t, t, u, MPFR_RNDN);  /* t = 2*pi*x/u * (1 + theta3)^3 where
+                                            |theta3| <= 2^-prec */
+      /* if t is zero here, it means the division by u underflowd */
+      if (MPFR_UNLIKELY (MPFR_IS_ZERO (t)))
+        {
+          mpfr_set_ui (y, 1, MPFR_RNDZ);
+          if (MPFR_IS_LIKE_RNDZ(rnd_mode,0))
+            {
+              inexact = -1;
+              mpfr_nextbelow (y);
+            }
+          else
+            inexact = 1;  
+          goto end;
+        }
+      /* since prec >= 2, |(1 + theta3)^3 - 1| <= 4*theta3 <= 2^(2-prec) */
+      expt = MPFR_GET_EXP (t);
+      /* we have |s| <= 2^(expt + 2 - prec) */
+      mpfr_cos (t, t, MPFR_RNDN);
+      err = expt + 2 - prec;
+      expt = MPFR_GET_EXP (t); /* new exponent of t */
+      /* the total error is at most 2^err + ulp(t)/2 = 2^err + 2^(expt-prec-1)
+         thus if err <= expt-prec-1, it is bounded by 2^(expt-prec),
+         otherwise it is bounded by 2^(err+1). */
+      err = (err <= expt - prec - 1) ? expt - prec : err + 1;
+      /* normalize err for mpfr_can_round */
+      err = expt - err;
+      if (MPFR_CAN_ROUND (t, err, precy, rnd_mode))
+        break;
+      /* check exact cases: this can only occur if 2*pi*x/u is a multiple
+         of pi/2, i.e., if x/u is a multiple of 1/4 */
+      if (nloops == 1)
+        {
+          inexact = mpfr_div_ui (t, x, u, MPFR_RNDZ);
+          mpfr_mul_2ui (t, t, 2, MPFR_RNDZ);
+          if (inexact == 0 && mpfr_integer_p (t))
+            {
+              if (mpfr_odd_p (t))
+                /* t is odd: we have kpi+pi/2, thus cosu = 0,
+                   for the sign, we always return +0, following IEEE 754-2019:
+                   cosPi(n + 1/2) is +0 for any integer n when n + 1/2 is
+                   representable. */
+                mpfr_set_zero (y, +1);
+              else /* t is even: case kpi */
+                {
+                  mpfr_div_2ui (t, t, 1, MPFR_RNDZ);
+                  if (!mpfr_odd_p (t))
+                    /* case 2kpi: cosu = 1 */
+                    mpfr_set_ui (y, 1, MPFR_RNDZ);
+                  else
+                    mpfr_set_si (y, -1, MPFR_RNDZ);
+                }
+              goto end;
+            }
+        }
+      MPFR_ZIV_NEXT (loop, prec);
+    }
+  MPFR_ZIV_FREE (loop);
+
+  inexact = mpfr_set (y, t, rnd_mode);
+
+ end:
+  mpfr_clear (t);
+  MPFR_SAVE_EXPO_FREE (expo);
+  return underflow ? inexact : mpfr_check_range (y, inexact, rnd_mode);
+}
diff --git a/src/mpfr.h b/src/mpfr.h
index 4e6959f61..74d6c675a 100644
--- a/src/mpfr.h
+++ b/src/mpfr.h
@@ -722,6 +722,7 @@ __MPFR_DECLSPEC int mpfr_csc (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
 __MPFR_DECLSPEC int mpfr_cot (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
 
 __MPFR_DECLSPEC int mpfr_sinu (mpfr_ptr, mpfr_srcptr, unsigned long, mpfr_rnd_t);
+__MPFR_DECLSPEC int mpfr_cosu (mpfr_ptr, mpfr_srcptr, unsigned long, mpfr_rnd_t);
 
 __MPFR_DECLSPEC int mpfr_hypot (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t);
 __MPFR_DECLSPEC int mpfr_erf (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 7d7f0e0b7..190caaf13 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -34,13 +34,13 @@ TESTS_NO_TVERSION = tabort_prec_max tassert tabort_defalloc1		\
      talloc-cache tasin tasinh tatan tatanh 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 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 tfmma tfmod    \
-     tfms tfpif tfprintf tfrac tfrexp tgamma tgamma_inc tget_d          \
-     tget_d_2exp tget_f tget_flt tget_ld_2exp tget_q tget_set_d64       \
-     tget_set_d128 tget_sj tget_str tget_z tgmpop tgrandom thyperbolic  \
-     thypot tinp_str							\
+     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 \
+     tfmma tfmod tfms tfpif tfprintf tfrac tfrexp tgamma tgamma_inc     \
+     tget_d tget_d_2exp tget_f tget_flt tget_ld_2exp tget_q             \
+     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	\
      tlog_ui tmin_prec tminmax tmodf tmul tmul_2exp tmul_d tmul_ui	\
      tnext tnrandom tnrandom_chisq tout_str toutimpl tpow tpow3		\
diff --git a/tests/tcosu.c b/tests/tcosu.c
new file mode 100644
index 000000000..5fa742f76
--- /dev/null
+++ b/tests/tcosu.c
@@ -0,0 +1,174 @@
+/* Test file for mpfr_cosu.
+
+Copyright 2020 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"
+
+static void
+test_singular (void)
+{
+  mpfr_t x, y;
+  int inexact;
+
+  mpfr_init (x);
+  mpfr_init (y);
+
+  /* check u = 0 */
+  mpfr_set_ui (x, 17, MPFR_RNDN);
+  inexact = mpfr_cosu (y, x, 0, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_nan_p (y));
+
+  /* check x = NaN */
+  mpfr_set_nan (x);
+  inexact = mpfr_cosu (y, x, 1, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_nan_p (y));
+
+  /* check x = +Inf */
+  mpfr_set_inf (x, 1);
+  inexact = mpfr_cosu (y, x, 1, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_nan_p (y));
+
+  /* check x = -Inf */
+  mpfr_set_inf (x, -1);
+  inexact = mpfr_cosu (y, x, 1, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_nan_p (y));
+
+  /* check x = +0 */
+  mpfr_set_zero (x, 1);
+  inexact = mpfr_cosu (y, x, 1, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0);
+  MPFR_ASSERTN(inexact == 0);
+
+  /* check x = -0 */
+  mpfr_set_zero (x, -1);
+  inexact = mpfr_cosu (y, x, 1, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0);
+  MPFR_ASSERTN(inexact == 0);
+
+  mpfr_clear (x);
+  mpfr_clear (y);
+}
+
+static void
+test_exact (void)
+{
+  mpfr_t x, y;
+  int inexact;
+
+  mpfr_init (x);
+  mpfr_init (y);
+
+  /* check 2*pi*x/u = pi/2 thus x/u = 1/4, for example x=1 and u=4 */
+  mpfr_set_ui (x, 1, MPFR_RNDN);
+  inexact = mpfr_cosu (y, x, 4, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_zero_p (y) && mpfr_signbit (y) == 0);
+  MPFR_ASSERTN(inexact == 0);
+
+  /* check 2*pi*x/u = pi thus x/u = 1/2, for example x=2 and u=4 */
+  mpfr_set_ui (x, 2, MPFR_RNDN);
+  inexact = mpfr_cosu (y, x, 4, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_cmp_si (y, -1) == 0 && inexact == 0);
+
+  /* check 2*pi*x/u = 3*pi/2 thus x/u = 3/4, for example x=3 and u=4 */
+  mpfr_set_ui (x, 3, MPFR_RNDN);
+  inexact = mpfr_cosu (y, x, 4, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_zero_p (y) && mpfr_signbit (y) == 0);
+  MPFR_ASSERTN(inexact == 0);
+
+  /* check 2*pi*x/u = 2*pi thus x/u = 1, for example x=4 and u=4 */
+  mpfr_set_ui (x, 4, MPFR_RNDN);
+  inexact = mpfr_cosu (y, x, 4, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0 && inexact == 0);
+
+  /* check 2*pi*x/u = -pi/2 thus x/u = -1/4, for example x=-1 and u=4 */
+  mpfr_set_si (x, -1, MPFR_RNDN);
+  inexact = mpfr_cosu (y, x, 4, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_zero_p (y) && mpfr_signbit (y) == 0);
+  MPFR_ASSERTN(inexact == 0);
+
+  /* check 2*pi*x/u = -pi thus x/u = -1/2, for example x=-2 and u=4 */
+  mpfr_set_si (x, -2, MPFR_RNDN);
+  inexact = mpfr_cosu (y, x, 4, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_cmp_si (y, -1) == 0 && inexact == 0);
+
+  /* check 2*pi*x/u = -3*pi/2 thus x/u = -3/4, for example x=-3 and u=4 */
+  mpfr_set_si (x, -3, MPFR_RNDN);
+  inexact = mpfr_cosu (y, x, 4, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_zero_p (y) && mpfr_signbit (y) == 0);
+  MPFR_ASSERTN(inexact == 0);
+
+  /* check 2*pi*x/u = -2*pi thus x/u = -1, for example x=-4 and u=4 */
+  mpfr_set_si (x, -4, MPFR_RNDN);
+  inexact = mpfr_cosu (y, x, 4, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0 && inexact == 0);
+
+  mpfr_clear (x);
+  mpfr_clear (y);
+}
+
+static void
+test_regular (void)
+{
+  mpfr_t x, y, z;
+  int inexact;
+
+  mpfr_init2 (x, 53);
+  mpfr_init2 (y, 53);
+  mpfr_init2 (z, 53);
+
+  mpfr_set_ui (x, 17, MPFR_RNDN);
+  inexact = mpfr_cosu (y, x, 42, MPFR_RNDN);
+  /* y should be cos(2*17*pi/42) rounded to nearest */
+  mpfr_set_str (z, "-0xd.38462625fd3ap-4", 16, MPFR_RNDN);
+  MPFR_ASSERTN(mpfr_equal_p (y, z));
+  MPFR_ASSERTN(inexact < 0);
+
+  mpfr_clear (x);
+  mpfr_clear (y);
+  mpfr_clear (z);
+}
+
+/* FIXME[VL]: For mpfr_cosu, the range reduction should not be expensive.
+   If I'm not mistaken, this is linear in the bitsize of the exponent
+   since one just needs to compute the argument modulo the integer u. */
+#define TEST_FUNCTION mpfr_cosu
+#define ULONG_ARG2
+#ifndef MPFR_USE_MINI_GMP
+#define REDUCE_EMAX 262143 /* otherwise arg. reduction is too expensive */
+#else
+#define REDUCE_EMAX 16383  /* reduce further since mini-gmp works in O(n^2) */
+#endif
+#include "tgeneric.c"
+
+int
+main (void)
+{
+  tests_start_mpfr ();
+
+  test_singular ();
+  test_exact ();
+  test_regular ();
+
+  test_generic (MPFR_PREC_MIN, 100, 1);
+
+  tests_end_mpfr ();
+  return 0;
+}