Finished rewrite of AGM in algorithms.tex, added bibliography entry.

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

2 files changed, 312 insertions, 135 deletions

diff --git a/doc/algorithms.bib b/doc/algorithms.bib
index b6ad6b3..b6af68e 100644
--- a/doc/algorithms.bib
+++ b/doc/algorithms.bib
@@ -1,3 +1,12 @@
+@article {Cox84,
+   author       = {David A. Cox},
+   title        = {The arithmetic-geometric mean of {G}auss},
+   journal      = {L'Enseignement Mathématique},
+   year         = {1984},
+   volume       = {30},
+   pages        = {275--330}
+}
+
 @phdthesis {Dupont06,
    author       = {R\'egis Dupont},
    title        = {Moyenne arithm\'etico-g\'eom\'etrique, suites de
diff --git a/doc/algorithms.tex b/doc/algorithms.tex
index 275ab3a..8b27b04 100644
--- a/doc/algorithms.tex
+++ b/doc/algorithms.tex
@@ -1,4 +1,4 @@
-\documentclass [12pt]{article}
+\documentclass [11pt]{article}
 
 \usepackage[a4paper]{geometry}
 \usepackage[utf8]{inputenc}
@@ -24,8 +24,8 @@
 \newcommand {\round}{\operatorname {\circ}}
 \DeclareMathOperator{\pinf}{\bigtriangleup}
 \DeclareMathOperator{\minf}{\bigtriangledown}
-\DeclareMathOperator{\N}{\mathcal N}
 \DeclareMathOperator{\A}{\mathcal A}
+\newcommand {\N}{\mathbb N}
 \newcommand {\Z}{\mathbb Z}
 \newcommand {\Q}{\mathbb Q}
 \newcommand {\R}{\mathbb R}
@@ -35,6 +35,7 @@
 \renewcommand {\leq}{\leqslant}
 \renewcommand {\geq}{\geqslant}
 \newcommand {\AGM}{\operatorname{AGM}}
+\newcommand {\sign}{\operatorname{sign}}
 
 \newtheorem{theorem}{Theorem}
 \newtheorem{lemma}[theorem]{Lemma}
@@ -50,7 +51,7 @@
 
 \title {MPC: Algorithms and Error Analysis}
 \author {Andreas Enge \and Philippe Th\'eveny \and Paul Zimmermann}
-\date {Draft; September 18, 2012}
+\date {Draft; January 23, 2014}
 
 \begin {document}
 \maketitle
@@ -538,6 +539,7 @@ For $0 < x_1 \leq 1$, we have $\appro {x_1}^2/2 \leq \appro {x_1}/2$ and
 \subsection {Complex functions}
 
 \subsubsection {Addition/subtraction}
+\label {sssec:propadd}
 
 Using the notation introduced above, we consider
 \[
@@ -1899,84 +1901,200 @@ Thus, assuming that $n 2^{-p} \leq 1$, the error estimates
 \eqref {eq:powui_re} and \eqref {eq:powui_im} are valid with $n$
 in the place of $n - 1$.
 
-\subsection{\texttt {mpc\_agm}}
+\subsection{\texttt {mpc\_agm}, \texttt {mpc\_agm1}}
+
+\paragraph {Definition.}
+Let $a$, $b$ be non-zero complex numbers.
+Define sequences of arithmetic means $(a_n)$
+and geometric means $(b_n)$
+by $\corr {a_0} = a$, $\corr {b_0} = b$,
+$\corr {a_{n+1}} = \frac {\corr {a_n} + \corr {b_n}}{2}$ and
+$\corr {b_{n+1}} = \sqrt {\corr {a_n} \corr {b_n}}$.
+At each step, there is a choice of sign for the square root.
+If $\corr {a_n}$ and $\corr {b_n}$ are at an angle different from $0$
+and $\pi$, then
+they define a two-dimensional pointed cone in the complex plane.
+Notice that $\corr {a_{n+1}}$ lies in this cone.
+Following \cite {Cox84} we call \emph {right} the choice that makes
+also $\corr {b_{n+1}}$ lie in the cone, and following \cite {CrTh10}
+we call the resulting sequences \emph {optimal}.
+There is a common limit of the sequences, the
+\emph {arithmetic-geometric mean}
+$\AGM (a, b)$.
+
+It is immediate that $\AGM$ is symmetric, that is,
+$\AGM (a, b) = \AGM (b, a)$, and homogeneous, that is,
+$\AGM (\lambda a, \lambda b) = \lambda \AGM (a, b)$ for any non-zero
+complex number $\lambda$.
+
+So we may assume that $|a| \geq |b|$, and
+$\AGM (a, b) = a \AGM (a_0, b_0)$
+with $\corr {a_0} = 1$, $\corr {b_0} = b/a$, $|\corr {b_0}| \leq 1$.
+
+We need to examine the corner cases.
+If one or both of $a$ and $b$ are zero, all geometric means are zero,
+and $\AGM (a, b) = 0$.
+If the angle between $a$ and $b$ is~$0$, then $\corr {b_0}$ is a positive real
+number, and $\AGM (1, \corr {b_0})$ may be computed with \mpfr.
+If the angle between $a$ and $b$ is $\pi$, then $\corr {b_0}$ is a negative real
+number larger than or equal to $-1$. The arithmetic mean of $1$ and $-1$
+is zero, so $\AGM (1, -1) = 0$.
+If $\corr {b_0} > -1$, we take the first geometric mean with a positive imaginary
+part, so that also $\Im (\AGM (1, \corr {b_0}))$ will be positive.
+Notice that $\Im (\AGM (1, \corr {b_0}))$ has the same sign as
+$\Im (\corr {b_0})$ unless
+$\corr {b_0}$ is real, so this choice determines our branch cut for $\AGM$.
+
+So in the following, we analyse the computation of $\AGM (1, \appro {b_0})$
+with $\appro {b_0} = \round (\corr {b_0}) \in \C \backslash \R^{\geq 0}$,
+where the real and imaginary parts of $\appro {b_0}$ are rounded towards~$0$,
+which ensures $| \appro {b_0} | \leq 1$
+with absolute errors of the real and imaginary parts at most \ulp {1}.
+In the following, relative errors will be most convenient to work with;
+by Proposition~\ref {prop:comabstorelerror} we have
+$\relerror (\appro {b_0}) \leq 2^{1-p}$, where $p$ is the working precision.
+
+
+\paragraph {The first iteration -- entering a quadrant.}
+If $\Re (\appro {b_0}) < 0$, then significant cancellation can occur for the
+arithmetic mean in the first iteration, which thus needs to be analysed
+separately.
+
+From now on, we use arbitrary rounding modes and apply
+Proposition~\ref {prop:comrelround} with $c = 1$
+(except for $\Re (\appro {a_1})$, which we need to round down or,
+equivalently, towards~$0$ for a later step in the analysis).
+We let $\corr {b_1} = \sqrt {\corr {b_0}}$ and
+$\appro {b_1} = \round \left( \sqrt {\appro {b_0}} \right)$ with
+\[
+\relerror (\appro {b_1}) \leq
+2^{1-p} + \left( 1 + 2^{1-p} \right) 2^{1-p}
+\leq 2^2 \cdot 2^{1-p}
+\]
+by \eqref {eq:propsqrt} and Proposition~\ref {prop:comrelround}.
+We let $\corr {a_1} = \frac {1 + \corr {b_0}}{2}$ and
+$\appro {a_1} = \round \left( \frac {1 + \appro {b_0}}{2} \right)$.
+The imaginary part of $\appro {a_1}$ has an error of at most \ulp{1},
+but following the discussion of \S\ref {sssec:propadd}, its real part
+may have an absolute error measured in $\Ulp$ of
+\[
+2^{\Exp (\Re (\appro {b_0})) - \Exp (\Re (\appro {a_1}))} + 1
+\leq
+2^{- \Exp (\Re (\appro {a_1}))} + 1
+\leq
+2^{- \Exp (\Re (\appro {a_1})) + 1},
+\]
+where we used that both $\Re (\appro {b_0})$ and $\Re (\appro {a_1})$ have
+non-positive exponents.
+By Proposition~\ref {prop:comabstorelerror},
+\[
+\relerror (\appro {a_1})
+\leq
+2^{- \Exp (\Re (\appro {a_1})) + 1} \cdot 2^{1-p}.
+\]
+If $\Re (\appro {b_0}) \geq 0$, then $\Exp (\Re (\appro {a_1})) = 0$, and
+we obtain a bound that is similar to the one above for
+$\relerror (\appro {b_1})$. If $\appro {b_0}$ approaches $-1$, then
+the relative error in $\appro {a_1}$ becomes arbitrarily bad, as
+$\appro {a_1}$ becomes arbitrarily small.
+
+Notice that, independently of the rounding mode,
+$\appro {a_1}$ and $\appro {b_1}$ lie in the same complex quadrant
+of numbers having non-negative real part and an imaginary part of the same
+sign as that of $\appro {b_0}$
+(or of positive imaginary part if $\appro {b_0}$ is real).
+During the remainder of the algorithm, we will not leave this quadrant
+and thus not see any more cancellation in the arithmetic mean.
+
+
+\paragraph {The idea of the analysis.}
+
+Let us first give the basic idea of the following, rather technical analysis.
+Assume a target precision of $N$ bits, that is, a target
+relative error of $2^{-N}$. If $a$ and $b$ are of roughly equal size, the AGM
+iteration converges quadratically, that is, the number of correct digits
+doubles in each step, and we need about $\log_2 N$ iterations. Unfortunately,
+when $a$ and $b$ are of very different sizes, we have slower convergence.
+To illustrate this, consider $a = 1$ and $b = 2^{-e}$; during the first
+iterations, $a_n \approx 2^{-n}$ and $b_n \approx 2^{-e/2^n}$. So we need
+to increase the number of iterations by roughly $\log_2$ of the absolute
+value of the exponent of $b/a$; the precise bound $B (N)$ on the number
+of iterations is given in~\eqref {eq:agmbound}.
+Moreover, when $\Re (b/a) < 0$, the situation of
+very different exponents in $a$ and $b$ may occur after one iteration through
+cancellation as explained above.
+
+Unlike Newton iterations, AGM iterations are not auto-correcting, but
+rounding errors propagate. We lose a constant number of bits per iteration;
+so to reach the desired precision of $N$ bits, we need to carry out all
+computations at a working precision of about $p = N + c B (N, b/a)$.
+The following discussion provides explicit bounds for all these quantities.
+
+
+\paragraph {Error propagation.}
 
 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 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})$.
-
-Among others, we will show in the following, using induction, that
-$\Re (\appro {a_n})$, $\Re (\appro {b_n}) \geq 0$ and
-$\Im (\appro {a_n})$, $\Im (\appro {b_n}) > 0$, regardless of the rounding mode.
-By assumption, this holds for $n = 0$.
-
-Let $c_n = \frac {\appro {a_{n-1}} + \appro {b_{n-1}}}{2}$, so that
-$\appro {a_n} = \round(c_n)$ satisfies
-$\Re (\appro {a_n}) \geq 0$, $\Im (\appro {a_n}) > 0$.
-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
+$\corr {a_n} = \frac {\corr {a_{n-1}} + \corr {b_{n-1}}}{2}$,
+$\corr {c_n} = \corr {a_{n-1}} \corr {b_{n-1}}$,
+$\corr {b_n} = \sqrt {\corr {c_n}}$.
+The computation of $\appro {a_1}$ and $\appro {b_1}$ and their error
+analysis have been given above.
+For $n \geq 2$, we compute the sequences
+\begin {eqnarray*}
+\appro {a_n}
+& = & \round \left( \frac {\appro {a_{n-1}} + \appro {b_{n-1}}}{2} \right), \\
+\appro {c_n}
+& = & \round \left( \appro {a_{n-1}} \appro {b_{n-1}} \right), \\
+\appro {b_n}
+& = & \round \left( \sqrt {\appro {c_n}} \right).
+\end {eqnarray*}
+Then one sees by induction that
+$\appro {a_n}$ and $\appro {b_n}$ lie in the same quadrant as
+$\appro {a_1}$ and $\appro {b_1}$
+(or, for that matter, $\corr {a_1}$ and $\corr {b_1}$).
+
+Let $\alpha_n = \relerror (\appro {a_n}) / 2^{1-p}$,
+$\gamma_n = \relerror (\appro {c_n}) / 2^{1-p}$ and
+$\beta_n = \relerror (\appro {b_n}) / 2^{1-p}$.
+
+By \eqref {eq:propaddrel} and Proposition~\ref {prop:comrelround},
 \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}
+\label {eq:agmalpha}
+\alpha_n \leq
+\sqrt 2 \, \max (\alpha_{n-1}, \beta_{n-1})
+\left( 1 + 2^{1-p} \right) + 1.
 \end {equation}
-by Proposition~\ref {prop:comrelround}.
-
-Let $d_n = \round \left( \appro {a_{n-1}} \appro {b_{n-1}} \right)$,
-so that $\Im (d_n) \geq 0$, and $\Im (d_n) = 0$ occurs only if
-$\appro {a_{n-1}}$ and $\appro {b_{n-1}}$ are both purely imaginary
-and $\Re (d_n) < 0$.
 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}.
+\label {eq:agmgamma}
+\gamma_n \leq
+\left( \alpha_{n-1} + \beta_{n-1} + \alpha_{n-1} \beta_{n-1} 2^{1-p} \right)
+\left( 1 + 2^{1-p} \right) + 1.
 \end {equation}
-Now $\appro {b_n} = \round (\sqrt {d_n})$ satisfies
-$\Re (\appro {b_n}) \geq 0$ and $\Im (\appro {b_n}) > 0$,
-and by \eqref {eq:propsqrt} and Proposition~\ref {prop:comrelround},
-assuming $\zeta_n \leq 1$, we obtain
+If we assume for the time being that $\gamma_n \leq 1 / (2^{1-p})$
+(which will be shown under additional constraints later), then
+\eqref {eq:propsqrt} and Proposition~\ref {prop:comrelround} imply
 \begin {equation}
-\label {eq:agmepsilon}
-\epsilon_n \leq
-\zeta_n + c (1 + \zeta_n) 2^{1-p}.
+\label {eq:agmbeta}
+\beta_n \leq
+\gamma_n (1 + 2^{1-p}) + 1.
 \end {equation}
 
-Let $r_n = (2^n - 1) d$ for some $d > 2 c$. Then one can show that
-for a sufficiently high
-precision~$p$ and not too many steps~$n$ compared to~$p$, one has
-$\eta_n, \epsilon_n \leq r_n 2^{1-p}$.
-Now we will prove this explicitly for $d=2^k$ with $k \geq 2$
-and $c \leq 1$, and show
-$\eta_n, \epsilon_n \leq r_n 2^{1-p} \leq 2^{n + k + 1 - p}$
-for $p \geq 2 (n + k) - 1$
-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
+Let $r_n = (2^n - 1) r_1$ for some integer $r_1$ such that
+$\alpha_1$, $\beta_1 \leq r_1$.
+We may use $r_1 = 2^k$ with $k = \max (2, - \Exp (\Re (\appro {a_1}) + 1))$
+by the discussion of the first iteration above.
+We now show by induction that
+$\alpha_n, \beta_n \leq r_n \leq 2^{n + k}$ for $p \geq 2 (n + k) - 1$,
+that is, if the number of iterations $n$ is not too large compared to
+the working precision $p$.
+By \eqref {eq:agmgamma} we obtain
+\[
+\gamma_n \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}.
 \]
-From $r_{n-1} < 2^{k + n-1}$ we obtain
+From $r_{n-1} < 2^{n - 1 + k}$ we deduce
 \[
 r_{n-1}^2 2^{1-p} < 2^{2 (n + k) - 1 - p} \leq 1
 \]
@@ -1984,33 +2102,35 @@ as long as $p \geq 2 (n + k) - 1$, and then, under the same condition,
 \[
 \left( r_{n-1}^2 + 2 r_{n-1} + r_{n-1}^2 2^{1-p} \right) 2^{1-p}
 \leq
-(r_{n-1} + 1)^2 \, 2^{1-p} \leq (2^{k+n-1})^2 \, 2^{1-p} \leq 1.
+(r_{n-1} + 1)^2 \, 2^{1-p} \leq (2^{n+k-1})^2 \, 2^{1-p} \leq 1.
 \]
 Then
 \[
-\frac {\zeta_n}{2^{1-p}}
+\gamma_n
 \leq
 2 r_{n-1} + 2
 \leq
 r_n - 2.
 \]
-By \eqref {eq:agmepsilon},
+By \eqref {eq:agmbeta},
 \[
-\frac {\epsilon_n}{2^{1-p}} \leq
+\beta_n
+\leq
 r_n - 1 + r_n 2^{1-p}
 \leq
 r_n - 1 + 2^{n+k+1-p}
 \leq r_n
 \]
-under a milder than the previous condition. Finally by \eqref {eq:agmeta},
+under a milder than the previous condition. Finally by \eqref {eq:agmalpha},
 for $p \geq 3$,
 \[
-\frac {\eta_n}{2^{1-p}} \leq
+\alpha_n
+\leq
 \sqrt 2 \, r_{n-1} + 1 + \sqrt {2} \, r_{n-1} 2^{1-p}
-\leq \frac{5}{4} \sqrt{2} r_{n-1} + 1 \leq 2 r_{n-1} + 2^k = r_n.
+\leq \frac{5}{4} \sqrt{2} r_{n-1} + 1 \leq 2 r_{n-1} + r_1 = r_n.
 \]
 This finishes the induction argument.
-To summarise, we have
+To summarise, we have for any fixed $n$ that
 \begin {equation}
 \label {eq:propagm}
 \corr {a_n} = (1 + \theta_1) \appro {a_n}
@@ -2020,40 +2140,123 @@ To summarise, we have
 2 (n + k) - 1 \leq p.
 \end {equation}
 
-We now use \cite[Prop.~3.3]{Dupont06}, which states that for
-$z_1 \neq 0, 1$ and
-\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
+
+\paragraph {Mathematical error.}
+
+We also need to estimate the error made by carrying out only a finite number
+of iterations. Let $z$ satisfy $\Re (z) \geq 0$ and $z \neq 0, 1$, and
+consider the optimal AGM sequences $(a_n')$ and $(b_n')$ computed (with
+infinite precision) from $a_0' = 1$ and $b_0' = z$. Let $N' \in \N$ and
+\[
+n \geq B' (N', z)
+  := \max \left( 1, \left\lceil \log_2 |\log_2 |z|| \right\rceil \right)
+        + \lceil \log_2 (N'+2) \rceil + 2
+\]
+(where $\log_2 0$ is to be understood as $- \infty$).
+Then by \cite[Prop.~3.3, p.~88]{Dupont06},
+\[
+a_n' = (1 + \theta_2) \AGM (1, z)
+\text { with }
+|\theta_2| \leq 2^{-N'}.
+\]
+(Notice that here, the relative error as defined in
+\cite[Def.~1.2, p.~20]{Dupont06} is taken with the roles of the correct
+and the approximated values reversed compared to our definition.)
+
+In our context, we may have $\Re (\corr {b_0}) < 0$, but after one iteration,
+$\Re (\corr {b_1} / \corr {a_1}) \geq 0$. So we consider
+$\corr {z} = b_0' = \corr {b_1} / \corr {a_1}$, so that by
+homogeneity, $\corr {a_n} = \corr {a_1} a_{n-1}'$,
+$\corr {b_n} = \corr {a_1} b_{n-1}'$ and
+$\AGM (1, \corr {b_0}) = \AGM (\corr {a_1}, \corr {b_1})
+= \corr {a_1} \AGM (1, \corr {z})$.
+Here
+\[
+|z|
+= \left| \frac {\corr {b_1}}{\corr {a_1}} \right|
+= \frac {|\sqrt {\corr {b_0}}|}{|\corr {a_1}|}
+= \frac {|\corr {b_0}|}{2 |\corr {a_1}|}
+\]
+is bounded above by
+\[
+\frac {1}{2 \, \Re (\corr {a_1})}
+\leq
+2^{- \Exp (\Re (\corr {a_1}))}
+=
+2^{- \Exp (\Re (\appro {a_1}))}
+\]
+(since $\Re (\appro {a_1})$ was rounded down)
+and below by
+\[
+\frac {|\corr {b_0}|}{2}
+\geq
+\frac {1}{2} \, \max \left( \Re (\corr {b_0}), |\Im (\corr {b_0})| \right)
+\geq
+\frac {1}{4} \,
+2^{\max \left( \Exp (\Re (\corr {b_0})), \Exp (\Im (\corr {b_0})) \right)}
+=
+2^{\max \left( \Exp (\Re (\appro {b_0})), \Exp (\Im (\appro {b_0})) \right)
+- 2}
+\]
+(since both parts of $\appro {b_0}$ were rounded down).
+So for $n \geq B (N)$ with
+\begin{equation}
+\label{eq:agmbound}
+\begin {aligned}
+B (N)
+  = & \left\lceil \log_2 \max (- \Exp (\Re (\appro {a_1})),
+                             - \Exp (\Re (\appro {b_0})) + 2,
+                             - \Exp (\Im (\appro {b_0})) + 2
+                             ) \right\rceil \\
+    & + \lceil \log_2 (N+4) \rceil + 3
+\end {aligned}
 \end{equation}
-(where $\log_2 0$ is to be understood as $- \infty$),
-we have
+(in which all exponents are non-positive)
+we reach
 \[
-a_n = (1 + \theta_2) z
+\corr {a_n} = (1 + \theta_2) \AGM (1, \corr {b_0})
 \text { with }
 |\theta_2| \leq 2^{-(N+2)}.
 \]
-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
+
+
+\paragraph {Total error and working precision}
+
+Combining with \eqref {eq:propagm} we obtain
+$\AGM (1, b_0) = \frac {1 + \theta_1}{1 + \theta_2} \appro {a_n}
+= (1 + \theta) \appro {a_n}$ with
 \[
 |\theta| \leq \frac {|\theta_1| + |\theta_2|}{|1 - |\theta_2||}
-\leq 2 \left( 2^{n+k+1-p} + 2^{-(N+2)} \right).
+\leq \frac {4}{3} \, \left( 2^{n+k+1-p} + 2^{-(N+2)} \right)
 \]
+for $N \geq 2$.
+So after $n = B (N)$ steps of the AGM iteration at a working precision
+of $p \geq N + n + k + 3$, we obtain $\appro {a_n}$ with a relative error
+bounded by $\epsilon = \frac {2}{3} \cdot 2^{-N}$.
+
+Finally, we let $\corr {z} = \AGM (a, b) = a \AGM (1, \corr {b_0})$ and
+$\appro {z} = \round (\corr {a} \appro {a_n})$, where $a$ is known
+exactly.
+By~\eqref {eq:propmulrel} and Proposition~\ref {prop:comabstorelerror},
+using that $N \geq 2$ implies $B (N) \geq 7$ and $p \geq 12$,
+this leads to a relative error bounded by $2^{-N}$.
 
-So after $n = B (N, z_1)$ steps of the AGM iteration at a working precision
-of $p = N + n + k + 3$, we obtain $\appro z = \appro {a_n}$ with a relative error
-bounded by $2^{-N}$.
 
-Note that with $p = N + n + k + 3$, the constraint $2 (n + k) + 1 \leq p$ means
-$n + k \leq N+2$. Depending on the value of $z_1$, one might have to take
-$N$ larger than the required accuracy to ensure \eqref{eq:agmbound} is
-fulfilled.
+\paragraph {Parameter choice.}
 
-Using Propositions~\ref {prop:comrelerror} and~\ref {prop:relerror}, this
+Recall that in our analysis, $k = \max (2, - \Exp (\Re (\appro {a_1} + 1)))$
+is given by the input data. We may then choose a target relative error~$N$,
+which determines the number of iterations~$n = B (N)$ by~\eqref {eq:agmbound}.
+Then the working precision may be chosen as
+\[
+p = \max (2 (n + k) - 1,
+          N + n + k + 3).
+\]
+
+Using Propositions~\ref {prop:comrelerror} and~\ref {prop:relerror}, the
 complex relative error may be translated into an error expressed in $\Ulp$.
-With $\appro {z} = \appro x + i \appro y$, let
+With $\appro {z} = \appro x + i \appro y$
+the computed approximation of $\corr {z} = \AGM (a, b)$, let
 $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
@@ -2064,43 +2267,8 @@ 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.
 So one should let $k' = \max (k_R, k_I)$, replace $N$ by $N + k'$
-and adapt the precision accordingly.
-
-If $\Im(z_1) < 0$, then we can use the fact that $\AGM(1,\overline z_1)$ is
-the complex
-conjugate of $\AGM(1,z_1)$, thus the same error analysis applies;
-and if $\Im(z_1) = 0$, we are computing a real AGM and can call the
-corresponding MPFR function.
-
-Now assume $z_1$ is the rounding of some complex number $z_0$ with a
-relative error at most $2^{1-p}$. Then we have to replace $\epsilon_0 = 0$
-by $\epsilon_0 = 2^{1-p}$ in the above proof.
-This gives
-\[ \zeta_1 \leq \epsilon_0 + c (1 + \epsilon_0) 2^{1-p}
-           \leq \left( 2 + 2^{1-p} \right) 2^{1-p} \leq \frac{9}{4} \, 2^{1-p}, \]
-and
-\[ \epsilon_1 \leq \zeta_1 + c (1 + \zeta_1) 2^{1-p}
-               \leq \left( \frac{9}{4} + 1 + \frac{9}{4} \, 2^{1-p} \right) 2^{1-p}
-               \leq 4 \cdot 2^{1-p} \]
-as long as $p \geq 3$. Thus the bound $\epsilon_1 \leq r_1 2^{1-p}$ still
-holds in that case.
-By \eqref {eq:agmeta}, also $\eta_1 \leq 4 \cdot 2^{1-p} = r_1 2^{1-p}$.
-
-\paragraph{The general case.}
-
-Assume now that $a, b$ are two arbitrary complex numbers, and we want to
-approximate $z = \AGM(a, b)$.
-We first have to precisely define $\AGM(a, b)$. At each step indeed there are
-two choices for $b_n = \sqrt{a_{n-1} b_{n-1}}$. We take a \emph{good} choice
-so that $|a_n - b_n| \leq |a_n + b_n|$; this corresponds to an \emph{optimal}
-AGM sequence \cite{CrTh10}. Unless $a_{n-1} / b_{n-1}$ is real, which happens
-exactly for $a / b$ real, this sequence is unique.
-It is shown in \cite{Dupont06} that the good choice is homogeneous with respect
-to multiplication of $a$ and $b$ by a non-zero complex number $\lambda$,
-so that $\AGM (\lambda a, \lambda b) = \lambda \AGM (a, b)$.
-In particular, $\AGM(a,b) = a \AGM(1,b/a) = b \AGM (1, a/b)$.
-We can thus restrict ourselves to the special case $\AGM(1, z)$ studied above
-with moreover $|z| \leq 1$.
+and adapt the precision and number of iterations accordingly.
+
 
 \bibliographystyle{acm}
 \bibliography{algorithms}