From 470d5dfe899d050c8f60015778983197de226ce8 Mon Sep 17 00:00:00 2001
From: Fredrik Johansson <fredrik.johansson@gmail.com>
Date: Thu, 19 Apr 2012 13:02:05 +0200
Subject: [PATCH] include error analysis

---
 arb/zeta_ui_euler_product.c | 19 ++++++++++++++++++-
 1 file changed, 18 insertions(+), 1 deletion(-)

diff --git a/arb/zeta_ui_euler_product.c b/arb/zeta_ui_euler_product.c
index 0731e1e7..9354a99b 100644
--- a/arb/zeta_ui_euler_product.c
+++ b/arb/zeta_ui_euler_product.c
@@ -24,7 +24,23 @@
 ******************************************************************************/
 
 #include "arb.h"
-#include "arb.h"
+
+/*
+Let P(a,b) = prod_{a <= p <= b} (1 - p^(-s)).
+Then 1/zeta(s) = P(a,M) * P(M+1,inf).
+
+According to the analysis in S. Fillebrown,
+"Faster Computation of Bernoulli Numbers", Journal of Algorithms 13,
+431-445 (1992), it holds for all s >= 6 and M >= 1 that
+(1/P(M+1,inf) - 1) <= 2 * M^(1-s) / (s/2 - 1).
+
+Writing 1/zeta(s) = P(a,M) * (1 - eps) and solving for eps gives
+
+1/(1-eps) <= 1 + 2 * M^(1-s) / (s/2 - 1), so we have
+eps <= 2 * M^(1-s) / (s/2 - 1).
+
+Since 0 < P(a,M) <= 1, this bounds the absolute error of 1/zeta(s).
+*/
 
 static __inline__ void
 arb_mul_ui(arb_t y, const arb_t x, ulong c)
@@ -68,6 +84,7 @@ arb_zeta_inv_ui_euler_product(arb_t z, ulong s)
     }
 
     prec = arb_prec(z);
+    /* heuristic */
     wp = prec + FLINT_BIT_COUNT(prec) + (prec/s) + 4;
 
     arb_init(t, wp);