From 2a7a4098219700fafd28a9e6ff9ed9a73dadfc7a Mon Sep 17 00:00:00 2001 From: zimmerma Date: Mon, 8 Oct 2012 13:40:30 +0000 Subject: [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 --- doc/algorithms.tex | 55 ++++++++++++++++++++++++++++++++++-------------------- 1 file 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)$. -- cgit v1.2.1