rootofunity: Reformulate part of the error analysis.

1 files changed, 22 insertions, 18 deletions

diff --git a/src/rootofunity.c b/src/rootofunity.c
index 50af1e3..0882126 100644
--- a/src/rootofunity.c
+++ b/src/rootofunity.c
@@ -164,29 +164,33 @@ mpc_rootofunity (mpc_ptr rop, unsigned long n, unsigned long k, mpc_rnd_t rnd)
    mpq_mul_2exp (kn, kn, 1); /* kn=2*k/n < 2 */
 
    do {
-      prec += mpc_ceil_log2 (prec) + 5;
+      prec += mpc_ceil_log2 (prec) + 5; /* prec >= 6 */
 
       mpfr_set_prec (t, prec);
       mpfr_set_prec (s, prec);
       mpfr_set_prec (c, prec);
 
-      mpfr_const_pi (t, MPFR_RNDN); /* error <= 0.5 ulp but since
-                                       ulp(t)=2^(2-prec), the absolute error
-                                       is bounded by 2^(1-prec), and the
-                                       relative error is bounded by
-                                       2^(1-prec)/pi <= 0.64*2^(-prec) */
-      mpfr_mul_q (t, t, kn, MPFR_RNDN); /* error <= 1.15 ulp(t) */
-      /* Indeed, the error is bounded by 0.64*2^(-prec)*pi*kn + 0.5 ulp [1].
-         Applying 2^(-prec)*|x| <= ulp(x) to x=pi*kn, we get a bound of:
-         0.64*ulp(pi*kn)+0.5 ulp.
-         Now since ulp(pi*kn) <= 2*ulp(t), we get: |t-pi*kn| <= 1.78*ulp(t).
-         If we plug this into [1] we get:
-         0.64*2^(-prec)*t + 0.64*2^(-prec)*1.78*ulp(t) + 0.5 ulp(t)
-         <= 0.64*ulp(t) + 0.64*2^(-prec)*1.78*ulp(t) + 0.5 ulp(t)
-         <= 1.15*ulp(t) for prec >= 7. */
+      mpfr_const_pi (t, MPFR_RNDN);
+         /* The absolute error is bounded by 0.5 ulp(t) = 2^(1-prec),
+            and with prec >= 6 we have 50/2^4 < t < 51/2^4,
+            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)
+               = |kn*t/u| * 1.14 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.5*2^{Exp (t) - Exp (s resp. c)} + 0.5) ulp
+         /* 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:
@@ -199,10 +203,10 @@ mpc_rootofunity (mpc_ptr rop, unsigned long n, unsigned long k, mpc_rnd_t rnd)
             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.5*2^66+0.5) ulp < 2^67 ulp.
+            (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.5*2^(3 - Exp (s or c)) + 0.5) ulp
+            (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. */
 
    }