rootsofunity: use mean value theorem for analysis (suggested by Damien Robert)

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/mpc/branches/rootsunity@1197 211d60ee-9f03-0410-a15a-8952a2c7a4e4

2 files changed, 13 insertions, 18 deletions

diff --git a/doc/algorithms.tex b/doc/algorithms.tex
index cee55af..790eca1 100644
--- a/doc/algorithms.tex
+++ b/doc/algorithms.tex
@@ -431,32 +431,27 @@ Let
 \[
 \appro x = \cos {\appro {x_1}}.
 \]
-Using the addition formula for $\cos$,
+By the mean value theorem, there is a $\xi$ between $x_1$ and $\appro {x_1}$
+such that
 \[
-\cos (a + b) = \cos (a) \cos (b) - \sin (a) \sin (b),
+\cos (x_1) - \cos (\appro {x_1}) = -\sin (\xi) (x_1 - \appro {x_1}),
 \]
-we obtain
-\begin {eqnarray*}
+so that
+\[
 \error (\appro x)
-& \leq & |\cos (x)| (1 - \cos (\error (\appro {x_1})))
-+ |\sin (x) \sin (\error (\appro {x_1}))| \\
-& \leq & 2 \error (\appro {x_1})
-\end {eqnarray*}
-since $|\sin (\delta)|$, $1 - \cos (\delta) \leq \delta$
-(one even has $1 - \cos (\delta) \leq \frac {1}{2} \delta^2$,
-but this does not fundamentally improve the error bound).
-
+\leq \error (\appro {x_1}).
+\]
 Taking the exponents into account, one obtains
 \begin {equation}
 \label {eq:proprealcos}
 \error (\appro x)
 \leq
-2 k \, 2^{\Exp (\appro {x_1}) - \Exp (\appro x)}
+k \, 2^{\Exp (\appro {x_1}) - \Exp (\appro x)}
 \, 2^{\Exp (\appro x) - p}.
 \end {equation}
 
 For the sine function, a completely analogous argument shows that
-\eqref {eq:proprealcos} still holds.
+\eqref {eq:proprealcos} also holds.
 
 
 
diff --git a/src/rootofunity.c b/src/rootofunity.c
index 8576f62..92555ee 100644
--- a/src/rootofunity.c
+++ b/src/rootofunity.c
@@ -67,12 +67,12 @@ mpc_rootofunity (mpc_ptr rop, unsigned long int n, mpc_rnd_t rnd)
       mpfr_mul_2ui (t, t, 1u, GMP_RNDN);
       mpfr_div_ui (t, t, n, GMP_RNDN); /* error 2*0.5+0.5=1.5 ulp */
       mpfr_sin_cos (s, c, t, GMP_RNDN);
-         /* error (3*2^{Exp (t) - Exp (s resp.c)} + 0.5) ulp
-            <= 12.5 ulp for n>=3                             */
+         /* error (1.5*2^{Exp (t) - Exp (s resp.c)} + 0.5) ulp
+            <= 6.5 ulp for n>=3                             */
    }
-   while (   !mpfr_can_round (c, prec - 4, GMP_RNDN, GMP_RNDZ,
+   while (   !mpfr_can_round (c, prec - 3, GMP_RNDN, GMP_RNDZ,
                  MPC_PREC_RE(rop) + (MPC_RND_RE(rnd) == GMP_RNDN))
-          || !mpfr_can_round (s, prec - 4, GMP_RNDN, GMP_RNDZ,
+          || !mpfr_can_round (s, prec - 3, GMP_RNDN, GMP_RNDZ,
                  MPC_PREC_IM(rop) + (MPC_RND_IM(rnd) == GMP_RNDN)));
 
    inex_re = mpfr_set (mpc_realref(rop), c, MPC_RND_RE(rnd));