[algorithms.tex] improved the error analysis for agm1

git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@1284 211d60ee-9f03-0410-a15a-8952a2c7a4e4

1 files changed, 35 insertions, 20 deletions

diff --git a/doc/algorithms.tex b/doc/algorithms.tex
index 3ab0849..1e928ad 100644
--- a/doc/algorithms.tex
+++ b/doc/algorithms.tex
@@ -244,7 +244,7 @@ it is easy to switch back and forth between absolute and relative errors.
 
 \begin {prop}
 \label {prop:relerror}
-Let $\appro x$ be representable at precision $p$.
+Let $\appro x$ be a real number representable at precision $p$.
 \begin {enumerate}
 \item
 If $\error (\appro x) \leq k \Ulp (\appro x)$,
@@ -593,7 +593,7 @@ One readily verifies that
 \]
 \end {proof}
 
-Assume now that $\corr {z_n} = (1 + \theta_n) \appro {z_n}$
+Assume now that $\corr {z_n} = (1 + \theta_n) \appro {z_n}$ for $n=1,2$
 with $\epsilon_n = |\theta_n|$ lie in the same quadrant.
 Then
 $\corr z = (1 + \theta) \appro z$
@@ -1032,7 +1032,6 @@ there is $0 < \xi < \epsilon_1$ with
 \leq \frac {1}{2 \sqrt {1 - \epsilon_1}} \, \epsilon_1.
 \end {equation}
 For instance $\epsilon \leq \epsilon_1$ for $\epsilon_1 \leq \frac {3}{4}$.
-
 We even have $\epsilon \leq \epsilon_1$ for $\epsilon_1 \leq 1$,
 since $1 - \sqrt{1-x}$ is bounded by $x$
 for $0 < x \leq 1$, which comes from $1 - x \leq \sqrt{1-x}$.
@@ -1936,7 +1935,7 @@ the error of $d_n$ (relative to $a_{n-1} b_{n-1}$) is bounded above by
 \end {equation}
 Now $\appro {b_n} = \round (\sqrt {d_n})$, and by \eqref {eq:propsqrt}
 and Proposition~\ref {prop:comrelround},
-assuming $\zeta_n \leq \frac {3}{4}$,
+assuming $\zeta_n \leq 1$,
 we obtain
 \begin {equation}
 \label {eq:agmepsilon}
@@ -1947,7 +1946,8 @@ we obtain
 Let $r_n = (2^n - 1) d$ for some $d > 2 c$. Then for a sufficiently high
 precision~$p$ and not too many steps~$n$ compared to~$p$, we have
 $\eta_n, \epsilon_n \leq r_n 2^{1-p}$.
-For instance, $c=1$ and $d=4$ is compatible with any rounding mode, and
+Now we will prove this explicitly in the case $c=1$ and $d=4$,
+which is enough to deal with all rounding modes, and
 we show
 $\eta_n, \epsilon_n \leq r_n 2^{1-p} \leq 2^{n + 3 -p}$
 by induction over \eqref {eq:agmeta}, \eqref {eq:agmzeta}
@@ -1958,11 +1958,15 @@ Inductively, by \eqref {eq:agmzeta} we obtain
 2 r_{n-1} + 1 +
 \left( r_{n-1}^2 + 2 r_{n-1} + r_{n-1}^2 2^{1-p} \right) 2^{1-p}
 \leq
-r_n - 3 + 2^{2n+4-p}
+r_n - 3 + 2^{2n+3-p}
 \leq
 r_n - 2
 \]
-as long as $2n+4 \leq p$. Then by \eqref {eq:agmepsilon},
+as long as $2n+3 \leq p$.\footnote{Indeed,
+$r_{n-1}^2 2^{1-p} \leq 2^{2n-2} \cdot 2^4 \cdot 2^{1-p} = 2^{2n+3-p} \leq 1$,
+thus $(r_{n-1}^2 + 2 r_{n-1} + r_{n-1}^2 2^{1-p}) 2^{1-p} \leq
+(r_{n-1} + 1)^2 2^{1-p} \leq 2^{2n-2} \cdot 2^4 \cdot 2^{1-p} = 2^{2n+3-p}$.}
+Then by \eqref {eq:agmepsilon},
 \[
 \frac {\epsilon_n}{2^{1-p}} \leq
 r_n - 1 + r_n 2^{1-p}
@@ -1971,12 +1975,12 @@ r_n - 1 + 2^{n+3-p}
 \leq r_n
 \]
 under a milder than the previous condition. Finally by \eqref {eq:agmeta},
+for $p \geq 3$,
 \[
 \frac {\eta_n}{2^{1-p}} \leq
 \sqrt 2 \, r_{n-1} + 1 + \sqrt {2} \, r_{n-1} 2^{1-p}
-\leq r_n
+\leq \frac{5}{4} \sqrt{2} r_{n-1} + 1 \leq 2 r_{n-1} + d = r_n.
 \]
-under an even less restrictive condition.
 
 To summarise, we have
 \begin {equation}
@@ -1990,11 +1994,11 @@ To summarise, we have
 
 We now use \cite[Prop.~3.3]{Dupont06}, which states that for
 $z_1 \neq 0, 1$,
-\[
+\begin{equation} \label{eq:agmbound}
 n \geq B (N, z_1)
   = \max \left( 1, \left\lceil \log_2 |\log_2 |z_1|| \right\rceil \right)
        + \lceil \log_2 (N+4) \rceil + 2
-\]
+\end{equation}
 (where $\log_2 0$ is to be understood as $- \infty$),
 we have
 \[
@@ -2015,26 +2019,37 @@ So after $n = B (N, z_1)$ steps of the AGM iteration at a working precision
 of $p = N + n + 5$, we obtain $\appro z = \appro {a_n}$ with a relative error
 bounded by $2^{-N}$.
 
+Note that with $p = N + n + 5$, the constraint $2n + 3 \leq p$ gives
+$n \leq N+2$. Depending on the value of $z_1$, one might have to take
+$N$ larger than the required accuracy to ensure Eq.~(\ref{eq:agmbound}) is
+fulfilled.
+
 Using Propositions~\ref {prop:comrelerror} and~\ref {prop:relerror}, this
 complex relative error may be translated into an error expressed in $\Ulp$.
 With $\appro {z} = \appro x + i \appro y$, let
-$k_R = \max (\Exp (\appro y) - \Exp (\appro x) + 1, 0) + 1$,
-$k_I = \max (\Exp (\appro x) - \Exp (\appro y) + 1, 0) + 1$ and
-$k = \max (k_R, k_I)$.
+$k_R = \max (\Exp (\appro y) - \Exp (\appro x) + 1, 0) + 1$, and
+$k_I = \max (\Exp (\appro x) - \Exp (\appro y) + 1, 0) + 1$.
 Then we have
 $\error (\appro x) \leq 2^{k_R - N + p} \Ulp (\appro x)$ and
 $\error (\appro y) \leq 2^{k_I - N + p} \Ulp (\appro y)$.
 
-In practice, one should take this additional loss into account: In a first
-loop, one may use an arbitrary guess for~$k$; if rounding fails after the
-first computation, one has nevertheless obtained the correct value of~$k$.
-Then after $n = B (N + k, z_1)$ steps of the AGM iteration at a working
+In practice, one should take this additional loss into account:
+if rounding fails after the
+first computation, nevertheless the values of $k_R$ and $k_I$ will most likely
+not change with a larger precision.
+Then take $k = \max (k_R, k_I)$,
+after $n = B (N + k, z_1)$ steps of the AGM iteration at a working
 precision of $p = N + k + n + 5$, one has
 $\error (\appro x)
-\leq 2^{p - N - (k - k_R)} \Ulp (\appro x) \leq 2^{p - N} \Ulp (\appro x)$
+\leq 2^{p - N - (k - k'_R)} \Ulp (\appro x)$
 and
 $\error (\appro y)
-\leq 2^{p - N - (k - k_I)} \Ulp (\appro y) \leq 2^{p - N} \Ulp (\appro y)$.
+\leq 2^{p - N - (k - k'_I)} \Ulp (\appro y)$,
+where $k'_R, k'_I$ denote the new values of $k_R$ and $k_I$.
+It then suffices to check a posteriori that $k'_R \leq k$ and
+$k'_I \leq k$, then
+$\error (\appro x) \leq 2^{p - N} \Ulp (\appro x)$ and
+$\error (\appro y) \leq 2^{p - N} \Ulp (\appro y)$.