rootofunity: Move the error analysis from the source code to the

documentation,

* src/rootofunity.c: Drop error analysis.
* doc/algorithms.tex" Add error analysis.

2 files changed, 76 insertions, 39 deletions

diff --git a/doc/algorithms.tex b/doc/algorithms.tex
index 128a1d9..a7b5bcc 100644
--- a/doc/algorithms.tex
+++ b/doc/algorithms.tex
@@ -51,7 +51,7 @@
 
 \title {MPC: Algorithms and Error Analysis}
 \author {Andreas Enge \and Philippe Th\'eveny \and Paul Zimmermann}
-\date {Draft; May 24, 2016}
+\date {Draft; November 15, 2017}
 
 \begin {document}
 \maketitle
@@ -1901,6 +1901,78 @@ 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\_rootofunity}}
+
+Consider co-prime integers $n > 0$ and $k \geq 0$. We wish to compute
+the primitive $n$-th root of unity $e^{2 \pi i k / n}$.
+First, we compute $\pi$ with rounding to nearest and an absolute error
+of \ulp{0.5}. Using Proposition~\ref {prop:relerror} we obtain a relative
+error of $\relerror (\appro {\pi}) \leq 2^{-p}$ at precision~$p$, which
+would be enough for the following analysis. However, we can be a bit more
+precise, knowing the first few binary digits of $\pi$: As soon as $p \geq 8$,
+we have $201/2^6 \leq \appro {\pi}$, so
+\[
+\relerror (\appro {\pi})
+= \left| \frac {\appro {\pi} - \corr {\pi}}{\appro {\pi}} \right|
+\leq 0.637 \cdot 2^{-p}.
+\]
+Let $r$ be the (exact) rational number $2k/n$.
+Next we compute
+$\appro {t} = \round (r \appro {\pi})$, as an approximation to
+$\corr {t} = r \pi$, rounded to nearest.
+Its absolute error is bounded as follows:
+\begin {align*}
+| \appro {t} - \corr {t} |
+& \leq |r \appro {\pi} - r \corr {\pi}| + 0.5 \, \Ulp (\appro {t}) \\
+& =  |r \appro {\pi} / \appro {t}| \cdot |\appro {t}|
+     \cdot |(\appro {\pi} - \corr {\pi}) / \appro {\pi}|
+     + 0.5 \, \Ulp (\appro {t}) \\
+& \leq |r \appro {\pi} / \appro {t}| \cdot 2^{\Exp (\appro {t})}
+     \cdot 0.637 \cdot 2^{-p} + 0.5 \, \Ulp (\appro {t}) \\
+& \leq \left( \max (|r \appro {\pi} / \appro {t}|, 1) \cdot 0.637
+              + 0.5 \right) \Ulp (\appro {t})
+  \text { by Definition~\ref {def:ulp}}.
+\end {align*}
+The factor $|r \appro {\pi} / \appro {t}|$ is larger than~$1$ only if
+rounding down has occurred for $|\appro {t}|$ by at most \ulp {0.5};
+with $p \geq 8$, the worst case bound is then $128.5/128$,
+and the absolute error is bounded by $1.15 \, \Ulp (\appro {t})$.
+
+The result is now obtained as
+\[
+\appro {s} + \appro {c} \, i
+\text { with }
+\appro {s} = \round (\sin (\appro {t}))
+\text { and }
+\appro {c} = \round (\cos (\appro {t})),
+\]
+rounded to nearest.
+By \eqref {eq:proprealcos}, the absolute error for the real part satisfies
+\begin {align*}
+|\appro {s} - \corr {s}|
+& \leq \left( 1.15 \cdot 2^{\Exp (\appro {t}) - \Exp (\appro {s})}
+   + 0.5 \right) \, \Ulp (\appro {s}) \\
+& \leq \left( 1.15 \cdot 2^{3 - \Exp (\appro {s})} + 0.5 \right)
+   \, \Ulp (\appro {s}) \text { since } \Exp (\appro {t}) \leq 3 \\
+& \leq 2^{4 - \Exp (\appro {s})} \, \Ulp (\appro {s})
+   \text { since } \Exp (\appro {s}) \leq 0,
+\end {align*}
+and the analogous bound holds for the imaginary part.
+
+If $n$ is large and $k$ close to a multiple of $n/4$, then $c$ or $s$
+can become small, and we lose about as many bits as given by the exponent
+of $s$ or $c$. The exact minimum of $s$ and~$c$ depends on the residue
+class of $n$ modulo~$4$. To simplify the argument, we consider the $4 n$-th
+roots of unity, which contain the $n$-th ones, and which are symmetrically
+distributed with respect to the real and imaginary axes. Then the smallest
+absolute value of the real or imaginary part is given by
+$\sin (2 \pi / (4 n))$. Using $\sin (x) > x - x^3 / 6$ for $x > 0$
+we obtain $\Exp (\appro {s}) \geq -63$ for $n < 2^{64}$; otherwise said,
+the number of bits lost during the computation is at most $67$ in this
+case, and in general at most $3$ more than the number of bits in~$n$.
+
+
 \subsection {\texttt {mpc\_cmp\_abs}}
 
 Let $z_1 = x_1 + i y_1$ and $z_2 = x_2 + i y_2$. We want to check whether
diff --git a/src/rootofunity.c b/src/rootofunity.c
index 9b4b7be..4b6974f 100644
--- a/src/rootofunity.c
+++ b/src/rootofunity.c
@@ -156,7 +156,9 @@ mpc_rootofunity (mpc_ptr rop, unsigned long n, unsigned long k, mpc_rnd_t rnd)
 
    prec = MPC_MAX_PREC(rop);
 
-   mpfr_init2 (t, 67); /* see the argument at the end of the following loop */
+   /* For the error analysis justifying the following algorithm,
+      see algorithms.tex. */
+   mpfr_init2 (t, 67);
    mpfr_init2 (s, 67);
    mpfr_init2 (c, 67);
    mpq_init (kn);
@@ -171,45 +173,8 @@ mpc_rootofunity (mpc_ptr rop, unsigned long n, unsigned long k, mpc_rnd_t rnd)
       mpfr_set_prec (c, prec);
 
       mpfr_const_pi (t, MPFR_RNDN);
-         /* The absolute error is bounded by 0.5 ulp(t) = 2^(1-prec),
-            and with prec >= 2 we have 50/2^4 <= t,
-            so the relative error is bounded above by
-            2^(1-prec)/t <= 0.64*2^(-prec); otherwise said,
-            |(t-pi)/t| <= 0.64*2^(-prec). */
       mpfr_mul_q (t, t, kn, MPFR_RNDN);
-         /* We analyse mpfr_mul_q (u, t, kn, MPFR_RNDN).
-            |u-kn*pi| <= |kn*t-kn*pi| + 0.5 ulp(u)
-               =  |kn*t/u| * |u| * |(t-pi)/t| + 0.5 ulp(u)
-               <= |kn*t/u| * 2^Exp(u) * 0.64*2^(-prec) + 0.5 ulp(u)
-               <= max(|kn*t/u|,1) * 1.14 ulp(u)
-                                           since 2^Exp(u)*2^(-prec) < ulp(u).
-               The factor |kn*t/u| is larger than 1 only if rounding down has
-               occurred for |u| by at most 0.5 ulp(u); with prec >= 6, the
-               worst case bound is then 32.5/32,
-               and the total error bounded by 1.16 ulp(u). */
-
       mpfr_sin_cos (s, c, t, MPFR_RNDN);
-         /* error (1.16*2^{Exp (t) - Exp (s resp. c)} + 0.5) ulp
-            according to Section 1.2.3 on error propagation of the sine
-            and cosine functions in algorithms.tex.
-            We have 0<t<2*pi, so Exp (t) <= 3.
-            Unfortunately s or c can be close to 0, but with n<2^64,
-            we lose at most about 64 bits:
-            Where the minimum of s and c over all primitive n-th roots of
-            unity is reached depends on n mod 4.
-            To simplify the argument, we will consider the 4*n-th roots of
-            unity, which contain the n-th roots of unity and which are
-            symmetrically distributed with respect to the real and imaginary
-            axes, so that it becomes enough to consider only s for k=1.
-            With n<2^64, the absolute values of all s or c are at least
-            sin (2 * pi * 2^(-64) / 4) > 2^(-64) of exponent at least -63.
-            So the error is bounded above by
-            (1.16*2^66+0.5) ulp < 2^67 ulp.
-            To obtain a more precise bound for smaller n, which is useful
-            especially at small precision, we may use the error bound of
-            (1.16*2^(3 - Exp (s or c)) + 0.5) ulp
-            < 2^(4 - Exp (s or c)) ulp, since Exp (s or c) is at most 0. */
-
    }
    while (   !mpfr_can_round (c, prec - (4 - mpfr_get_exp (c)),
                  MPFR_RNDN, MPFR_RNDZ,