Merges changesets 4629 and 4630 from the trunk (bug fixes, updated

algorithms.tex).


git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/2.3@4631 280ebfd0-de03-0410-8827-d642c229c3f4

3 files changed, 23 insertions, 33 deletions

diff --git a/algorithms.tex b/algorithms.tex
index ce14d412a..bf0f702d4 100644
--- a/algorithms.tex
+++ b/algorithms.tex
@@ -870,7 +870,9 @@ above reasoning. At step (1c) we have $r=y/2$, thus $r-y$ simplifies to $-r$.
 \subsection{The cosine function}
 
 To evaluate $\cos x$ with a target precision of $n$ bits, we use the following
-algorithm with working precision $m$:
+algorithm with working precision $m$, after an additive argument reduction
+which reduces $x$ in the interval $[-\pi, \pi]$, using the
+\texttt{mpfr\_remainder} function:
 \begin{quote}
 $k \leftarrow \lfloor \sqrt{n/2} \rfloor$ \\
 $r \leftarrow x^2$ rounded up \\ % err <= ulp(r)
@@ -907,7 +909,7 @@ giving $\epsilon_k \le (2 l_0+1/3) 2^{2k-m}$.
 \subsection{The sine function}
 
 The sine function is computed from the cosine, with a working precision of
-$m$ bits:
+$m$ bits, after an additive argument reduction in $[-\pi, \pi]$:
 \begin{quote}
 $c \leftarrow \cos x$ rounded away \\
 $t \leftarrow c^2$ rounded away \\
@@ -993,34 +995,20 @@ bound of $2^{-p+7}$.
 
 \subsection{The tangent function}
 
-The tangent function is computed from the cosine, 
-using $\tan x = {\rm sign}(x) \sqrt{\frac{1}{\cos^2 x} - 1}$,
+The tangent function is computed from the \texttt{mpfr\_sin\_cos} function,
+which computes simultaneously $\sin x$ and $\cos x$
 with a working precision of $m$ bits:
 \begin{quote}
-$c \leftarrow \cos x$ rounded down \\ % c <= cos(x) <= 1
-$t \leftarrow c^2$ rounded down \\    % t <= cos(x)^2 <= 1
-$v \leftarrow 1/t$ rounded up \\      % v >= 1/cos(x)^2 >= 1
-$u \leftarrow v - 1$ rounded up \\    % u >= 1/cos(x)^2 - 1
-$s \leftarrow {\rm sign}(x) \sqrt{u}$ rounded away from $0$ \\
+$s, c \leftarrow \circ(\sin x), \circ(\cos x)$ \quad [to nearest] \\
+$t \leftarrow \circ(s/c)$ \quad [to nearest] \\
 \end{quote}
-The absolute error on $c$ is at most $\ulp(c)$.
-Hence the error on
-$t$ is at most $\ulp(t) + 2 c \ulp(c) \le 5 \ulp(t)$,
-% err(t) <= ulp(t) + |c^2-c'^2| <= ulp(t) + |c-c'|*|c+c'|
-%        <= ulp(t) + 2*ulp(c)*c <= 5*ulp(t)
-that on $v$ is at most $\ulp(v) + 5 \ulp(t)/t^2
-        \le \ulp(v) + 10 \ulp(1/t) \le 11 \ulp(v)$,
-% err(v) <= ulp(v) + |1/t-1/t'| <= ulp(v) + |t-t'|/t/t'
-%        <= ulp(v) + 5*ulp(t)/t^2 <= ulp(v) + 10*ulp(1/t) <= 11*ulp(v)
-that on $u$ is at most $\ulp(u) + 11 \ulp(v) \le (1 + 11 \cdot 2^e) \ulp(u)$
-where $e$ is the exponent difference between $v$ and $u$.
-% err(u) <= ulp(u) + err(v) <= ulp(v) + 11*ulp(v) <= (1+11*2^e)*ulp(u)
-The final error on $s$ is $\le \ulp(s) + (1+11 \cdot 2^e) \ulp(u)/2/\sqrt{u}
-        \le \ulp(s) + (1+11 \cdot 2^e) \ulp(u/\sqrt{u})
-        \le (2 + 11 \cdot 2^e) \ulp(s)$.
-% err(s) <= ulp(s) + |u-u'|/2/sqrt(u) <= ulp(s) + (1+11*2^e)*ulp(u)/2/sqrt(u)
-%        <= ulp(s) + (1+11*2^e)*ulp(u/sqrt(u))
-%        <= (2+11*2^e)*ulp(s) 
+We have $s = \sin(x) (1 + \theta_1)$ and $c = \cos(x) (1 + \theta_2)$
+with $|\theta_1|, |\theta_2| \leq 2^{-m}$, thus
+$t = (\tan x) (1 + \theta)^3$ with $|\theta| \leq 2^{-m}$.
+For $m \geq 2$, $|\theta| \leq 1/4$,
+$|(1 + \theta)^3 - 1| \leq 4 |\theta|$, thus we can write
+$t = (\tan x) (1 + 4 \theta)$, thus
+$|t - \tan x| \leq 4 \ulp(t)$.
 
 \subsection{The exponential function}
 
diff --git a/cos.c b/cos.c
index 35efff794..5ddcac8a7 100644
--- a/cos.c
+++ b/cos.c
@@ -231,10 +231,11 @@ mpfr_cos (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
         {
           if (m > k && (m - k >= precy + (rnd_mode == GMP_RNDN)))
             {
-              /* if round to nearest or away, result is s,
-                 otherwise it is round(nexttoward (s, 0)) */
-              if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (s)))
-                mpfr_nexttozero (s);
+              /* If round to nearest or away, result is s = 1 or -1,
+                 otherwise it is round(nexttoward (s, 0)). However in order to
+                 have the inexact flag correctly set below, we set |s| to
+                 1 - 2^(-m) in all cases. */
+              mpfr_nexttozero (s);
               break;
             }
         }
diff --git a/tan.c b/tan.c
index 4fe9b3a50..d38e5c1c9 100644
--- a/tan.c
+++ b/tan.c
@@ -62,6 +62,7 @@ mpfr_tan (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
   /* Compute initial precision */
   precy = MPFR_PREC (y);
   m = precy + MPFR_INT_CEIL_LOG2 (precy) + 13;
+  MPFR_ASSERTD (m >= 2); /* needed for the error analysis in algorithms.tex */
 
   MPFR_GROUP_INIT_2 (group, m, s, c);
   MPFR_ZIV_INIT (loop, m);
@@ -70,9 +71,9 @@ mpfr_tan (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
       /* The only way to get an overflow is to get ~ Pi/2
          But the result will be ~ 2^Prec(y). */
       mpfr_sin_cos (s, c, x, GMP_RNDN); /* err <= 1/2 ulp on s and c */
-      mpfr_div (c, s, c, GMP_RNDN);     /* err <= 2 ulps */
+      mpfr_div (c, s, c, GMP_RNDN);     /* err <= 4 ulps */
       MPFR_ASSERTD (!MPFR_IS_SINGULAR (c));
-      if (MPFR_LIKELY (MPFR_CAN_ROUND (c, m-1, precy, rnd_mode)))
+      if (MPFR_LIKELY (MPFR_CAN_ROUND (c, m - 2, precy, rnd_mode)))
         break;
       MPFR_ZIV_NEXT (loop, m);
       MPFR_GROUP_REPREC_2 (group, m, s, c);