Synchronize with r1286 from trunk.

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

4 files changed, 429 insertions, 53 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 6ca7454..669d651 100644
--- a/doc/algorithms.tex
+++ b/doc/algorithms.tex
@@ -3,6 +3,7 @@
 \usepackage[a4paper]{geometry}
 \usepackage[utf8]{inputenc}
 \usepackage[T1]{fontenc}
+\usepackage{ae}
 \usepackage{amsmath,amssymb}
 \usepackage{hyperref}
 \usepackage{comment}
@@ -31,6 +32,7 @@
 \renewcommand {\theta}{\vartheta}
 \renewcommand {\leq}{\leqslant}
 \renewcommand {\geq}{\geqslant}
+\newcommand {\AGM}{\operatorname{AGM}}
 
 \newtheorem{theorem}{Theorem}
 \newtheorem{lemma}[theorem]{Lemma}
@@ -46,7 +48,7 @@
 
 \title {MPC: Algorithms and Error Analysis}
 \author {Andreas Enge \and Philippe Th\'eveny \and Paul Zimmermann}
-\date {Draft; June 16, 2010}
+\date {Draft; September 18, 2012}
 
 \begin {document}
 \maketitle
@@ -242,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)$,
@@ -338,6 +340,50 @@ For the converse direction, write
 \]
 \end {proof}
 
+As a corollary to Propositions~\ref {prop:relerror}
+and~\ref {prop:comrelerror}, we obtain the following result that shows how
+to transform absolute into complex relative errors.
+
+\begin {prop}
+\label {prop:comabstorelerror}
+Let $\appro z = \appro x + i \appro y$ be representable at precision~$p$.
+If
+$\error (\appro x) \leq k \, \Ulp (\appro x)$ and
+$\error (\appro y) \leq k \, \Ulp (\appro y)$, then
+\[
+\relerror (\appro z) \leq k \, 2^{1-p}.
+\]
+As in Proposition~\ref {prop:relerror}, there is an immediate generalisation
+if $\appro z$ is not representable at precision~$p$.
+\end {prop}
+
+This result can be used to analyse how rounding affects complex
+relative errors.
+
+\begin {prop}
+\label {prop:comrelround}
+Let $\relerror (\appro z) = \epsilon$.
+Then
+\[
+\relerror (\round (\appro z)) \leq
+\epsilon + c (1 + \epsilon) 2^{1-p},
+\]
+where $c = \frac {1}{2}$ if both the real and the imaginary part are
+rounded to nearest, and $c = 1$ otherwise.
+\end {prop}
+
+\begin {proof}
+Write $\corr z = (1 + \theta) \appro z$ with $|\theta| = \epsilon$.
+Applying Proposition~\ref {prop:comabstorelerror} with $\appro z$ in the
+place of $\corr z$, $\round (\appro z)$ in the place of $\appro z$
+and $c$ in the place of~$k$,
+there is a $\theta'$ with $\appro z = (1 + \theta') \round (\appro z)$
+and $|\theta'| \leq c \, 2^{1-p}$.
+So $\corr z = (1 + \theta'') \round (\appro z)$ with
+$\theta'' = (1 + \theta)(1 + \theta') - 1
+= \theta + \theta' (1 + \theta)$.
+\end {proof}
+
 
 \subsection {Real functions}
 
@@ -424,6 +470,68 @@ Then by \cite[\S1.7]{MPFRAlgorithms},
 \end {equation}
 
 
+\subsubsection {Cosine and sine}
+\label {sssec:proprealcossin}
+
+Let
+\[
+\appro x = \cos {\appro {x_1}}.
+\]
+By the mean value theorem, there is a $\xi$ between $x_1$ and $\appro {x_1}$
+such that
+\[
+\cos (x_1) - \cos (\appro {x_1}) = -\sin (\xi) (x_1 - \appro {x_1}),
+\]
+so that
+\[
+\error (\appro x)
+\leq \error (\appro {x_1}).
+\]
+Taking the exponents into account, one obtains
+\begin {equation}
+\label {eq:proprealcos}
+\error (\appro x)
+\leq
+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} also holds.
+
+
+\subsubsection {Logarithm}
+\label {sssec:propreallog}
+
+Let
+\[
+\appro x = \log (1 + \appro {x_1})
+\]
+for $\appro {x_1} > -1$.
+By the mean value theorem, there is a $\xi$ between $x_1$ and $\appro {x_1}$
+such that
+\[
+\error (\appro x) = \frac {1}{1 + \xi} \error (\appro {x_1})
+\leq \frac {1}{1 + \min (x_1, \appro {x_1})} \error (\appro {x_1}).
+\]
+For $x_1 > 0$, this implies
+\begin {eqnarray*}
+\error (\appro x)
+& \leq & \error (\appro {x_1})
+\leq
+k \, 2^{\Exp (\appro {x_1}) - \Exp (\appro x)}
+\, 2^{\Exp (\appro x) - p} \\
+& \leq & 2 \, k \, \frac {\appro {x_1}}{\appro x} \, 2^{\Exp (\appro x) - p} \\
+& \leq & 2 \, k \, \frac {\appro {x_1}}{\appro {x_1} - \appro {x_1}^2/2}
+\, 2^{\Exp (\appro x) - p}
+\end {eqnarray*}
+using $\log (1 + z) \geq z - z^2/2$ for $z > 0$.
+For $0 < x_1 \leq 1$, we have $\appro {x_1}^2/2 \leq \appro {x_1}/2$ and
+\[
+\error (\appro x)
+\leq 4 \, k \, 2^{\Exp (\appro x) - p}.
+\]
+
 
 \subsection {Complex functions}
 
@@ -454,7 +562,7 @@ exponent than the operands, that is, if cancellation occurs.
 If $\appro {x_1}$ and $\appro {x_2}$, have the same sign, then there
 is no cancellation, $d_{R, n} \leq 0$ and
 \[
-\error (\appro x) \leq (k_{R,1} + k_{R,2} + c_R) 2^{\Exp (\appro x) - p}.
+\error (\appro x) \leq (k_{R,1} + k_{R,2}) 2^{\Exp (\appro x) - p}.
 \]
 
 An analogous error bound holds for the imaginary part.
@@ -463,6 +571,51 @@ For subtraction, the same bounds are obtained, except that the simpler bound
 now holds whenever $\appro {x_1}$ and $\appro {x_2}$ resp.
 $\appro {y_1}$ and $\appro {y_2}$ have different signs.
 
+We obtain a useful result for complex relative errors when $z_1$ and $z_2$
+lie in the same quadrant, so that cancellation occurs neither for the real
+nor for the imaginary part.
+
+\begin {lemma}
+\label {lm:arithgeom}
+Let $c_1 = a_1 + i b_1$ and $c_2 = a_2 + i b_2$ lie in the same quadrant,
+that is, $a_1 a_2$, $b_1 b_2 \geq 0$. Then
+\[
+|c_1| + |c_2| \leq \sqrt 2 \cdot |c_1 + c_2|.
+\]
+\end {lemma}
+
+\begin {proof}
+One readily verifies that
+\[
+2 |c_1 + c_2|^2 - (|c_1| + |c_2|)^2
+= (|c_1| - |c_2|)^2 + 4 (a_1 a_2 + b_1 b_2)
+\geq 0
+\]
+\end {proof}
+
+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$
+with
+\[
+\theta = \frac {\theta_1 \appro {z_1} + \theta_2 \appro {z_2}}
+               {\appro {z_1} + \appro {z_2}}.
+\]
+and
+\begin {equation}
+\label {eq:propaddrel}
+\relerror (\appro z)
+\leq
+\max (\epsilon_1, \epsilon_2)
+   \frac {|\appro {z_1}| + |\appro {z_2}|}{|\appro {z_1} + \appro {z_2}|}
+\leq
+\sqrt 2 \, \max (\epsilon_1, \epsilon_2)
+\end {equation}
+by Lemma~\ref {lm:arithgeom}.
+
+
+
 
 \subsubsection {Multiplication}
 \label {sssec:propmul}
@@ -605,7 +758,11 @@ A coarser bound may be obtained more easily by considering complex
 relative errors. Write $\corr {z_n} = (1 + \theta_n) \appro {z_n}$
 with $\epsilon_n = | \theta_n |$. Then $\corr z = (1 + \theta) \appro z$
 with $\theta = \theta_1 + \theta_2 + \theta_1 \theta_2$ and
-$\epsilon = |\theta| \leq \epsilon_1 + \epsilon_2 + \epsilon_1 \epsilon_2$.
+\begin {equation}
+\label {eq:propmulrel}
+\epsilon = \relerror (\appro z)
+\leq \epsilon_1 + \epsilon_2 + \epsilon_1 \epsilon_2.
+\end {equation}
 By Proposition~\ref {prop:relerror},
 we have $\epsilon_{X, n} \leq k_{X, n} 2^{1-p}$ for $X \in \{ R, I \}$,
 and by Proposition~\ref {prop:comrelerror},
@@ -844,6 +1001,42 @@ we find the exact same error estimate \eqref {eq:propmulcomrel}
 also for the case of division.
 
 
+\subsubsection {Square root}
+Let
+\[
+\appro z = \sqrt {\appro {z_1}}.
+\]
+Write $\corr {z_1} = (1 + \theta_1) \appro {z_1}$ with
+$\epsilon_1 = |\theta_1|$, and assume $\epsilon_1 < 1$.
+Then $\corr z = (1 + \theta) \appro z$ with
+\[
+\theta = \sqrt {1 + \theta_1} - 1
+= \frac {1}{2} \theta_1
++ \sum_{n=2}^\infty \frac {(-1)^{n+1} 1 \cdot 3 \cdots (2 n - 3)}{2^n \, n!}
+   \theta_1^n
+\]
+as a Taylor series, and
+\[
+\epsilon = |\theta|
+\leq
+\frac {1}{2}  \epsilon_1
++ \sum_{n=2}^\infty \frac {1 \cdot 3 \cdots (2 n - 3)}{2^n \, n!}
+\epsilon_1^n
+= 1 - \sqrt {1 - \epsilon_1}.
+\]
+By the mean value theorem, applied to the function $f (x) = \sqrt {1 - x}$,
+there is $0 < \xi < \epsilon_1$ with
+\begin {equation}
+\label {eq:propsqrt}
+\epsilon = \frac {1}{2 \sqrt {1 - \xi}} \, \epsilon_1
+\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}$.
+
+
 \subsubsection {Logarithm}
 
 
@@ -1699,7 +1892,201 @@ 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\_agm1}}
 
+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})$.
+
+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 1$,
+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}$.
+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}
+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+3-p}
+\leq
+r_n - 2
+\]
+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}
+\leq
+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 \frac{5}{4} \sqrt{2} r_{n-1} + 1 \leq 2 r_{n-1} + d = r_n.
+\]
+
+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}
+\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$,
+\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
+\[
+a_n = (1 + \theta_2) z
+\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
+\[
+|\theta| \leq \frac {|\theta_1| + |\theta_2|}{|1 - |\theta_2||}
+\leq 2 \left( 2^{n+3-p} + 2^{-(N+2)} \right)
+\]
+(see the computation at the end of \S\ref {sssec:propdiv}).
+
+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$, 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:
+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)$
+and
+$\error (\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)$.
+
+If $\Im(z_1) < 0$, then we can use the fact that $\AGM(1,\bar{z_1})$ is the
+conjugate of $\AGM(1,z_1)$, thus the same error analysis applies;
+and if $\Im(z_1) = 0$, we are computing a real AGM, we 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 (2 + 2^{1-p}) 2^{1-p} \leq \frac{9}{4} 2^{1-p}, \]
+and
+\[ \epsilon_1 \leq \zeta_1 + c (1 + \zeta_1) 2^{1-p}
+               \leq (\frac{9}{4} + 1 + \frac{9}{4} 2^{1-p}) 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.
+
+\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}.
+In case $|a_n - b_n| = |a_n + b_n|$, we choose $b_n$ such that
+$b_n/a_n$ has a positive imaginary part.
+
+In \cite{Dupont06}, it is shown that if $\lambda$ is a non-zero complex
+number, then $(\lambda a, \lambda b)$ is a good choice if and only if
+$(a, b)$ is a good choice. As a consequence, taking $\lambda = 1/a$,
+it shows that $\AGM(a,b) = a \AGM(1,b/a)$. We can thus restrict ourselves to
+the special case $\AGM(1, z)$ studied above.
+Since $\AGM(a,b) = \AGM(b,a)$, without loss of generality we can assume
+$|z| \leq 1$.
 
 \bibliographystyle{acm}
 \bibliography{algorithms}
diff --git a/doc/mpc.texi b/doc/mpc.texi
index a724877..33f1b03 100644
--- a/doc/mpc.texi
+++ b/doc/mpc.texi
@@ -5,7 +5,7 @@
 @synindex tp fn
 
 @set MINGMP 4.3.2
-@set MINMPFR 2.4.2
+@set MINMPFR 3.0.0
 
 @set AUTHORS Andreas Enge, Philippe Th@'eveny, Paul Zimmermann
 
@@ -13,7 +13,7 @@
 This manual is for GNU MPC, a library for multiple precision complex arithmetic,
 version @value{VERSION} of @value{UPDATED-MONTH}.
 
-Copyright @copyright{} 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 INRIA
+Copyright @copyright{} 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 INRIA
 
 @quotation
 Permission is granted to copy, distribute and/or modify this document
@@ -858,6 +858,11 @@ Set @var{rop} to @var{op1} @minus{} @var{op2} rounded according to @var{rnd}.
 For @code{mpc_ui_ui_sub}, @var{op1} is @var{re1} + @var{im1}.
 @end deftypefn
 
+@deftypefun int mpc_neg (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd})
+Set @var{rop} to @minus{}@var{op} rounded according to @var{rnd}.
+Just changes the sign if @var{rop} and @var{op} are the same variable.
+@end deftypefun
+
 @deftypefun int mpc_mul (mpc_t @var{rop}, mpc_t @var{op1}, mpc_t @var{op2}, mpc_rnd_t @var{rnd})
 @deftypefunx int mpc_mul_ui (mpc_t @var{rop}, mpc_t @var{op1}, unsigned long int @var{op2}, mpc_rnd_t @var{rnd})
 @deftypefunx int mpc_mul_si (mpc_t @var{rop}, mpc_t @var{op1}, long int @var{op2}, mpc_rnd_t @var{rnd})
@@ -877,6 +882,11 @@ in both cases rounded according to @var{rnd}.
 Set @var{rop} to the square of @var{op} rounded according to @var{rnd}.
 @end deftypefun
 
+@deftypefun int mpc_fma (mpc_t @var{rop}, mpc_t @var{op1}, mpc_t @var{op2}, mpc_t @var{op3}, mpc_rnd_t @var{rnd})
+Set @var{rop} to @var{op1}*@var{op2}+@var{op3},
+rounded according to @var{rnd}, with only one final rounding.
+@end deftypefun
+
 @deftypefun int mpc_div (mpc_t @var{rop}, mpc_t @var{op1}, mpc_t @var{op2}, mpc_rnd_t @var{rnd})
 @deftypefunx int mpc_div_ui (mpc_t @var{rop}, mpc_t @var{op1}, unsigned long int @var{op2}, mpc_rnd_t @var{rnd})
 @deftypefunx int mpc_div_fr (mpc_t @var{rop}, mpc_t @var{op1}, mpfr_t @var{op2}, mpc_rnd_t @var{rnd})
@@ -885,11 +895,6 @@ Set @var{rop} to the square of @var{op} rounded according to @var{rnd}.
 Set @var{rop} to @var{op1}/@var{op2} rounded according to @var{rnd}.
 @end deftypefun
 
-@deftypefun int mpc_neg (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd})
-Set @var{rop} to @minus{}@var{op} rounded according to @var{rnd}.
-Just changes the sign if @var{rop} and @var{op} are the same variable.
-@end deftypefun
-
 @deftypefun int mpc_conj (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd})
 Set @var{rop} to the conjugate of @var{op} rounded according to @var{rnd}.
 Just changes the sign of the imaginary part
@@ -907,24 +912,22 @@ Set the floating-point number @var{rop} to the norm of @var{op}
 rounded in the direction @var{rnd}.
 @end deftypefun
 
-@deftypefun int mpc_mul_2exp (mpc_t @var{rop}, mpc_t @var{op1}, unsigned long int @var{op2}, mpc_rnd_t @var{rnd})
+@deftypefun int mpc_mul_2ui (mpc_t @var{rop}, mpc_t @var{op1}, unsigned long int @var{op2}, mpc_rnd_t @var{rnd})
+@deftypefunx int mpc_mul_2si (mpc_t @var{rop}, mpc_t @var{op1}, long int @var{op2}, mpc_rnd_t @var{rnd})
 Set @var{rop} to @var{op1} times 2 raised to @var{op2}
-rounded according to @var{rnd}. Just increases the exponents
+rounded according to @var{rnd}. Just modifies the exponents
 of the real and imaginary parts by @var{op2}
 when @var{rop} and @var{op1} are identical.
 @end deftypefun
 
-@deftypefun int mpc_div_2exp (mpc_t @var{rop}, mpc_t @var{op1}, unsigned long int @var{op2}, mpc_rnd_t @var{rnd})
+@deftypefun int mpc_div_2ui (mpc_t @var{rop}, mpc_t @var{op1}, unsigned long int @var{op2}, mpc_rnd_t @var{rnd})
+@deftypefunx int mpc_div_2si (mpc_t @var{rop}, mpc_t @var{op1}, long int @var{op2}, mpc_rnd_t @var{rnd})
 Set @var{rop} to @var{op1} divided by 2 raised to @var{op2}
-rounded according to @var{rnd}. Just decreases the exponents
+rounded according to @var{rnd}. Just modifies the exponents
 of the real and imaginary parts by @var{op2}
 when @var{rop} and @var{op1} are identical.
 @end deftypefun
 
-@deftypefun int mpc_fma (mpc_t @var{rop}, mpc_t @var{op1}, mpc_t @var{op2}, mpc_t @var{op3}, mpc_rnd_t @var{rnd})
-Set @var{rop} to @var{op1} @times @var{op2} plus @var{op3},
-rounded according to @var{rnd}, with only one final rounding.
-@end deftypefun
 
 @node Power Functions and Logarithm
 @section Power Functions and Logarithm
diff --git a/doc/version.texi b/doc/version.texi
index b314fbd..7f5dbd4 100644
--- a/doc/version.texi
+++ b/doc/version.texi
@@ -1,4 +1,4 @@
-@set UPDATED 28 June 2012
-@set UPDATED-MONTH June 2012
-@set EDITION 1.0.0dev
-@set VERSION 1.0.0dev
+@set UPDATED 29 October 2012
+@set UPDATED-MONTH October 2012
+@set EDITION 1.1dev
+@set VERSION 1.1dev