algorithms.tex: AGM continued

algorithms.bib: references added and removed


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

2 files changed, 129 insertions, 57 deletions

diff --git a/doc/algorithms.bib b/doc/algorithms.bib
index 50bcaeb..2232a0e 100644
--- a/doc/algorithms.bib
+++ b/doc/algorithms.bib
@@ -1,3 +1,13 @@
+@phdthesis {Dupont06,
+   author       = {R\'egis Dupont},
+   title        = {Moyenne arithm\'etico-g\'eom\'etrique, suites de
+                   {B}orchardt et applications},
+   school       = {\'Ecole polytechnique},
+   address      = {Palaiseau},
+   type         = {Th\`ese de doctorat},
+   year         = {2006}
+}
+
 @Article{Friedland67,
   author = 	 {Paul Friedland},
   title = 	 {Algorithm 312: {A}bsolute value and square root of a complex
@@ -10,6 +20,12 @@
   annote = "\url{http://portal.acm.org/citation.cfm?id=363780}"
 }
 
+@Unpublished{MPFRAlgorithms,
+  author =       {{MPFR Team}},
+  title =        {The {MPFR} {L}ibrary: {A}lgorithms {A}nd {P}roofs},
+  note =         {\url{http://www.mpfr.org/algorithms.pdf}},
+}
+
 @Article{Smith98,
   author = 	 {David M. Smith},
   title = 	 {Algorithm 786: {M}ultiple-Precision Complex Arithmetic and
@@ -20,33 +36,3 @@
   number =	 4,
   pages =	 {359--367}
 }
-
-@Book{BrDiJeLeMeMuReStTo09,
-  author =       {Nicolas Brisebarre and Florent de Dinechin and
-                  Claude-Pierre Jeannerod and Vincent Lef\`evre and
-                  Guillaume Melquiond and Jean-Michel Muller and Nathalie
-                  Revol and Damien Stehl\'e and Serge Torres},
-  title =        {Handbook of Floating-Point Arithmetic},
-  publisher =    {Birkh\"auser},
-  year =         2009,
-  note =         {571 pages}
-}
-
-@Unpublished{MPFRAlgorithms,
-  author =       {{MPFR Team}},
-  title =        {The {MPFR} {L}ibrary: {A}lgorithms {A}nd {P}roofs},
-  note =         {\url{http://www.mpfr.org/algorithms.pdf}},
-}
-
-@Article{Percival03,
-  author = 	 "Colin Percival",
-  title = 	 "Rapid multiplication modulo the sum and difference of
-		  highly composite numbers",
-  journal =	 MC,
-  year =	 2003,
-  volume =	 72,
-  number =	 241,
-  pages =	 "387--395",
-  comment =	 "Includes good error analysis for FFT but proof is
-		  incorrect (result is OK)"
-}
diff --git a/doc/algorithms.tex b/doc/algorithms.tex
index 32f2254..90df53b 100644
--- a/doc/algorithms.tex
+++ b/doc/algorithms.tex
@@ -1891,34 +1891,120 @@ in the place of $n - 1$.
 
 \subsection{The AGM}
 
-We assume we compute the univariate AGM defined by $\AGM(1,z)$.
-It is easy to see
-that after one AGM iteration the values
-$\round \left( \frac {1+z}{2} \right)$ and
-$\round (\sqrt{z})$ are in the same quadrant (in the first quadrant if the
-imaginary part of $z$ is nonnegative, in the fourth quadrant otherwise).
-Taking into account symmetries, we thus assume $z$ is in the first quadrant.
-
-Let $\corr {a_n}, \corr {b_n}$ be the values of the AGM iteration performed
-with infinite precision, and
+Let
+\[
+z = \AGM (1, z_1);
+\]
+for the time being, we assume $\Re (z_1) \geq 0$ and $\Im (z_1) > 0$.
+
+Let $\corr {a_0} = \appro {a_0} = 1$,
+$\corr {b_0} = \appro {b_0} = z_1$, and let
+$\corr {a_n} = \frac {\corr {a_{n-1}} + \corr {b_{n-1}}}{2}$ and
+$\corr {b_n} = \sqrt {\corr {a_{n-1}} \corr {b_{n-1}}}$
+be the values of the AGM iteration performed with infinite precision and
 $\appro {a_n}, \appro {b_n}$ those performed with $p$-bit precision
-and rounding towards $+\infty$. We have $\corr {a_0} = a_0 = 1$ and
-$\corr {b_0} = b_0 = z$ (assuming $z$ is exactly representable in precision
-$p$).
-
-Each iteration computes
-$\appro {a_n} = \round \left(
-\frac {\appro {a_{n-1}} + \appro {b_{n-1}}}{2} \right)$ and
-$\appro {b_n} = \round \left(
-\sqrt {\appro {a_{n-1}} \appro {b_{n-1}}} \right)$.
-Assume by induction we can write
-$\appro {a_{n-1}} = \corr {a_{n-1}} (1 + \mu_{n-1})^{e_{n-1}}$
-and $\appro {b_{n-1}} = \corr {b_{n-1}} (1 + \nu_{n-1})^{e_{n-1}}$
-with $|\mu_{n-1}|, |\nu_{n-1}| < 2^{1-p}$.
-Then using the above lemma, since the division by $2$ is exact, we can write
-$\appro {a_n} = \corr {a_n} (1 + \mu_n)^{e_{n-1} + 1}$, and
-$\appro {b_n} = \corr {b_n} (1 + \nu_n)^{e_{n-1} + 2}$.
-Thus $e_n \leq e_{n-1} + 2$, which gives $e_n \leq 2n$.
+and some fixed rounding mode; let $c = \frac {1}{2}$ if this rounding mode
+is to nearest, $c = 1$ otherwise.
+Let $\eta_n = \relerror (\appro {a_n})$ and
+$\epsilon_n = \relerror (\appro {b_n})$.
+
+Let $c_n = \frac {\appro {a_{n-1}} + \appro {b_{n-1}}}{2}$, so that
+$\appro {a_n} = \round(c_n)$. By \eqref {eq:propaddrel}, the error of
+$c_n$ (relative to the correct value $a_n$) is bounded above by
+$\sqrt 2 \, \max (\epsilon_{n-1}, \eta_{n-1})$, division by~$2$ being exact.
+After rounding, we obtain
+\begin {equation}
+\label {eq:agmeta}
+\eta_n \leq
+\sqrt 2 \, \max (\eta_{n-1}, \epsilon_{n-1})
++ c \left( 1 + \sqrt 2 \, \max (\eta_{n-1}, \epsilon_{n-1}) \right) 2^{1-p}
+\end {equation}
+by Proposition~\ref {prop:comrelround}.
+
+Let $d_n = \round \left( \appro {a_{n-1}} \appro {b_{n-1}} \right)$.
+Then by \eqref {eq:propmulrel} and Proposition~\ref {prop:comrelround},
+the error of $d_n$ (relative to $a_{n-1} b_{n-1}$) is bounded above by
+\begin {equation}
+\label {eq:agmzeta}
+\zeta_n =
+\eta_{n-1} + \epsilon_{n-1} + \eta_{n-1} \epsilon_{n-1}
++ c (1 + \eta_{n-1} + \epsilon_{n-1} + \eta_{n-1} \epsilon_{n-1}) 2^{1-p}.
+\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}$,
+we obtain
+\begin {equation}
+\label {eq:agmepsilon}
+\epsilon_n \leq
+\zeta_n + c (1 + \zeta_n) 2^{1-p}.
+\end {equation}
+
+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
+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}
+and~\eqref {eq:agmepsilon}. For $n = 0$, we have $\eta_0 = \epsilon_0 = 0$.
+Inductively, by \eqref {eq:agmzeta} we obtain
+\[
+\frac {\zeta_n}{2^{1-p}} \leq
+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}
+\leq
+r_n - 2
+\]
+as long as $2n+4 \leq p$. Then by \eqref {eq:agmepsilon},
+\[
+\frac {\epsilon_n}{2^{1-p}} \leq
+r_n - 1 + r_n 2^{1-p}
+\leq
+r_n - 1 + 2^{n+3-p}
+\leq r_n
+\]
+under a milder than the previous condition. Finally by \eqref {eq:agmeta},
+\[
+\frac {\eta_n}{2^{1-p}} \leq
+\sqrt 2 \, r_{n-1} + 1 + \sqrt {2} \, r_{n-1} 2^{1-p}
+\leq r_n
+\]
+under an even less restrictive condition.
+
+To summarise, we have
+\begin {equation}
+\label {eq:propagm}
+\corr {a_n} = (1 + \theta_1) \appro {a_n}
+\text { with }
+|\theta_1| \leq 2^{n + 3 - p} \leq 2^{1 - p/2}
+\text { for }
+2n+4 \leq p.
+\end {equation}
+
+We now use \cite[Prop.~3.3]{Dupont06}, which states that for
+$z_1 \neq 0, 1$,
+\[
+n \geq B (N, z_1)
+  = \max \left( 1, \left\lceil \log_2 |\log_2 |z_1|| \right\rceil \right)
+       + \lceil \log_2 (N+2) \rceil + 2
+\]
+(where $\log_2 0$ is to be understood as $- \infty$),
+we have
+\[
+a_n = (1 + \theta_2) z
+\text { with }
+|\theta_2| \leq 2^{-N}.
+\]
+So taking $\appro z = \appro {a_n}$ as an approximation for $z$, we obtain
+$z = \frac {1 + \theta_1}{1 + \theta_2} \appro z
+= (1 + \theta) \appro z$ with
+\[
+|\theta| \leq \frac {|\theta_1| + |\theta_2|}{|1 - |\theta_2||}
+\]
+(see the computation at the end of \S\ref {sssec:propdiv}).
 
 
 \bibliographystyle{acm}