[sin.c,cos.c] use mpfr_sincos_fast when prec >= MPFR_SINCOS_THRESHOLD

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

6 files changed, 53 insertions, 13 deletions

diff --git a/cos.c b/cos.c
index cd4d6acd1..b0aeaa939 100644
--- a/cos.c
+++ b/cos.c
@@ -23,6 +23,16 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
 #define MPFR_NEED_LONGLONG_H
 #include "mpfr-impl.h"
 
+static int
+mpfr_cos_fast (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
+{
+  int inex;
+
+  inex = mpfr_sincos_fast (NULL, y, x, rnd_mode);
+  inex = inex >> 2; /* 0: exact, 1: rounded up, 2: rounded down */
+  return (inex == 2) ? -1 : inex;
+}
+
 /* f <- 1 - r/2! + r^2/4! + ... + (-1)^l r^l/(2l)! + ...
    Assumes |r| < 1/2, and f, r have the same precision.
    Returns e such that the error on f is bounded by 2^e ulps.
@@ -154,6 +164,13 @@ mpfr_cos (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
 
   /* Compute initial precision */
   precy = MPFR_PREC (y);
+
+  if (precy >= MPFR_SINCOS_THRESHOLD)
+    {
+      MPFR_SAVE_EXPO_FREE (expo);
+      return mpfr_cos_fast (y, x, rnd_mode);
+    }
+
   K0 = __gmpfr_isqrt (precy / 3);
   m = precy + 2 * MPFR_INT_CEIL_LOG2 (precy) + 2 * K0;
 
diff --git a/mpfr-impl.h b/mpfr-impl.h
index 7a83462f9..842b63069 100644
--- a/mpfr-impl.h
+++ b/mpfr-impl.h
@@ -1686,6 +1686,8 @@ __MPFR_DECLSPEC void mpfr_rand_raw _MPFR_PROTO((mp_ptr, gmp_randstate_t,
 __MPFR_DECLSPEC mpz_t* mpfr_bernoulli_internal _MPFR_PROTO((mpz_t*,
                                                             unsigned long));
 
+__MPFR_DECLSPEC int mpfr_sincos_fast _MPFR_PROTO((mpfr_t, mpfr_t, mpfr_srcptr, mp_rnd_t));
+
 #if defined (__cplusplus)
 }
 #endif
diff --git a/mpfr.texi b/mpfr.texi
index fcf518306..a25c797d6 100644
--- a/mpfr.texi
+++ b/mpfr.texi
@@ -1687,10 +1687,9 @@ cotangent of @var{op}, rounded in the direction @var{rnd}.
 @end deftypefun
 
 @deftypefun int mpfr_sin_cos (mpfr_t @var{sop}, mpfr_t @var{cop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
-Set simultaneously @var{sop} to the sine of @var{op} and
-                   @var{cop} to the cosine of @var{op},
-rounded in the direction @var{rnd} with the corresponding precisions of
-@var{sop} and @var{cop}, which must be different variables.
+Set simultaneously @var{sop} to the sine of @var{op} and @var{cop} to the
+cosine of @var{op}, rounded in the direction @var{rnd} with the corresponding
+precisions of @var{sop} and @var{cop}, which must be different variables.
 Return 0 iff both results are exact, more precisely it returns @math{s+4c}
 where @math{s=0} if @var{sop} is exact, @math{s=1} if @var{sop} is larger
 than the sine of @var{op}, @math{s=2} if @var{sop} is smaller than the sine
diff --git a/sin.c b/sin.c
index 37b953cbd..8fb9899a7 100644
--- a/sin.c
+++ b/sin.c
@@ -23,6 +23,16 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
 #define MPFR_NEED_LONGLONG_H
 #include "mpfr-impl.h"
 
+static int
+mpfr_sin_fast (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
+{
+  int inex;
+
+  inex = mpfr_sincos_fast (y, NULL, x, rnd_mode);
+  inex = inex & 3; /* 0: exact, 1: rounded up, 2: rounded down */
+  return (inex == 2) ? -1 : inex;
+}
+
 int
 mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
 {
@@ -62,6 +72,10 @@ mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
 
   /* Compute initial precision */
   precy = MPFR_PREC (y);
+
+  if (precy >= MPFR_SINCOS_THRESHOLD)
+    return mpfr_sin_fast (y, x, rnd_mode);
+
   m = precy + MPFR_INT_CEIL_LOG2 (precy) + 13;
   expx = MPFR_GET_EXP (x);
 
diff --git a/sin_cos.c b/sin_cos.c
index 707a30fa8..3842ad9d5 100644
--- a/sin_cos.c
+++ b/sin_cos.c
@@ -26,8 +26,6 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
 #define INEXPOS(y) ((y) == 0 ? 0 : (((y) > 0) ? 1 : 2))
 #define INEX(y,z) (INEXPOS(y) | (INEXPOS(z) << 2))
 
-int mpfr_sincos_fast (mpfr_ptr, mpfr_ptr, mpfr_srcptr, mp_rnd_t);
-
 /* (y, z) <- (sin(x), cos(x)), return value is 0 iff both results are exact
    ie, iff x = 0 */
 int
@@ -555,7 +553,11 @@ sincos_aux (mpfr_t s, mpfr_t c, mpfr_srcptr x, mp_rnd_t rnd_mode)
   return err;
 }
 
-/* assumes x is neither NaN, +/-Inf, nor +/- 0 */
+/* Assumes x is neither NaN, +/-Inf, nor +/- 0.
+   One of s and c might be NULL, in which case the corresponding value is
+   not computed.
+   Assumes s differs from c.
+ */
 int
 mpfr_sincos_fast (mpfr_t s, mpfr_t c, mpfr_srcptr x, mp_rnd_t rnd)
 {
@@ -565,7 +567,13 @@ mpfr_sincos_fast (mpfr_t s, mpfr_t c, mpfr_srcptr x, mp_rnd_t rnd)
   mp_exp_t err, errs, errc;
   MPFR_ZIV_DECL (loop);
 
-  w = MPFR_PREC(s) >= MPFR_PREC(c) ? MPFR_PREC(s) : MPFR_PREC(c);
+  MPFR_ASSERTN(s != c);
+  if (s == NULL)
+    w = MPFR_PREC(c);
+  else if (c == NULL)
+    w = MPFR_PREC(s);
+  else
+    w = MPFR_PREC(s) >= MPFR_PREC(c) ? MPFR_PREC(s) : MPFR_PREC(c);
   w += MPFR_INT_CEIL_LOG2(w) + 9; /* ensures w >= 10 (needed by sincos_aux) */
   mpfr_init2 (ts, w);
   mpfr_init2 (tc, w);
@@ -630,15 +638,15 @@ mpfr_sincos_fast (mpfr_t s, mpfr_t c, mpfr_srcptr x, mp_rnd_t rnd)
       /* adjust errors with respect to absolute values */
       errs = err - MPFR_EXP(ts);
       errc = err - MPFR_EXP(tc);
-      if (MPFR_CAN_ROUND (ts, w - errs, MPFR_PREC(s), rnd) &&
-          MPFR_CAN_ROUND (tc, w - errc, MPFR_PREC(c), rnd))
+      if ((s == NULL || MPFR_CAN_ROUND (ts, w - errs, MPFR_PREC(s), rnd)) &&
+          (c == NULL || MPFR_CAN_ROUND (tc, w - errc, MPFR_PREC(c), rnd)))
         break;
       MPFR_ZIV_NEXT (loop, w);
     }
   MPFR_ZIV_FREE (loop);
 
-  inexs = mpfr_set (s, ts, rnd);
-  inexc = mpfr_set (c, tc, rnd);
+  inexs = (s == NULL) ? 0 : mpfr_set (s, ts, rnd);
+  inexc = (c == NULL) ? 0 : mpfr_set (c, tc, rnd);
 
   mpfr_clear (ts);
   mpfr_clear (tc);
diff --git a/tuneup.c b/tuneup.c
index b000095fc..da221c50d 100644
--- a/tuneup.c
+++ b/tuneup.c
@@ -147,7 +147,7 @@ int verbose;
 
 
 /* First we include all the functions we want to tune inside this program.
-   We can't use GNU MPFR library since the THRESHOLD can't vary */
+   We can't use the GNU MPFR library since the thresholds are fixed macros. */
 
 /* Setup mpfr_exp_2 */
 mp_prec_t mpfr_exp_2_threshold;