Merge pull request flintlib#2138 from fredrik-johansson/zeta

Allow evaluating deflated Riemann zeta function and derivatives on small intervals containing 1

src/acb_poly/test/t-zeta_cpx_series.c

@@ -83,5 +83,62 @@ TEST_FUNCTION_START(acb_poly_zeta_cpx_series, state)
         _acb_vec_clear(z2, len);
     }
 
+    /* check deflations around 1 */
+    for (iter = 0; iter < 100 * 0.1 * flint_test_multiplier(); iter++)
+    {
+        acb_t s1, s2, a;
+        acb_ptr z1, z2;
+        slong i, len, prec;
+        int deflate;
+
+        acb_init(s1);
+        acb_init(s2);
+        acb_init(a);
+
+        acb_one(a);
+
+        prec = 2 + n_randint(state, 300);
+        len = 1 + n_randint(state, 20);
+        deflate = 1;
+
+        acb_randtest(s1, state, 300, 2);
+        acb_mul_2exp_si(s1, s1, -(slong) n_randint(state, 100));
+
+        acb_zero(s2);
+        arb_get_mag(arb_radref(acb_realref(s2)), acb_realref(s1));
+        arb_get_mag(arb_radref(acb_imagref(s2)), acb_imagref(s1));
+
+        acb_add_ui(s1, s1, 1, prec);
+        acb_add_ui(s2, s2, 1, prec);
+
+        z1 = _acb_vec_init(len);
+        z2 = _acb_vec_init(len);
+
+        _acb_poly_zeta_cpx_series(z1, s1, a, deflate, len, prec);
+        _acb_poly_zeta_cpx_series(z2, s2, a, deflate, len, prec);
+
+        for (i = 0; i < len; i++)
+        {
+            if (!acb_overlaps(z1 + i, z2 + i))
+            {
+                flint_printf("FAIL: overlap (deflation)\n\n");
+                flint_printf("iter = %wd\n", iter);
+                flint_printf("deflate = %d, len = %wd, i = %wd\n\n", deflate, len, i);
+                flint_printf("s1 = "); acb_printd(s1, prec / 3.33); flint_printf("\n\n");
+                flint_printf("s2 = "); acb_printd(s2, prec / 3.33); flint_printf("\n\n");
+                flint_printf("a = "); acb_printd(a, prec / 3.33); flint_printf("\n\n");
+                flint_printf("z1 = "); acb_printd(z1 + i, prec / 3.33); flint_printf("\n\n");
+                flint_printf("z2 = "); acb_printd(z2 + i, prec / 3.33); flint_printf("\n\n");
+                flint_abort();
+            }
+        }
+
+        acb_clear(a);
+        acb_clear(s1);
+        acb_clear(s2);
+        _acb_vec_clear(z1, len);
+        _acb_vec_clear(z2, len);
+    }
+
     TEST_FUNCTION_END(state);
 }

src/acb_poly/zeta_series.c

@@ -39,6 +39,66 @@ _acb_poly_zeta_cpx_series(acb_ptr z, const acb_t s, const acb_t a, int deflate,
         return;
     }
 
+    /* F(s) = zeta(s) - 1/(s-1) on small intervals containing 1 */
+    if (deflate && acb_is_one(a) && !acb_is_one(s) &&
+        (arb_contains_zero(acb_imagref(s)) && arb_contains_si(acb_realref(s), 1)))
+    {
+        acb_t t;
+        mag_t r, u, v;
+        slong n;
+
+        /* An extremely crude bound is |gamma_n| <= 8^(-n) n!.
+
+            (F^{(n)}(1+x) - F^{(n)}(1)) / n! = sum_{k>=1} F^{(n+k)} / (k! n!)
+            |(F^{(n)}(1+x) - F^{(n)}(1)) / n!| <= sum_{k>=1} x^k |gamma_{n+k}| / (k! n!)
+                                               <= sum_{k>=1} x^k (n+k)! / (k! n!) / 8^(n+k)
+                                               <= sum_{k>=1} x^k (n+1)^k / 8^(n+k)
+                                                = 8^(-n) (n+1) x / (8 - (n+1) x)
+        */
+
+        mag_init(r);
+        mag_init(u);
+        mag_init(v);
+        acb_init(t);
+
+        is_real = acb_is_real(s);
+
+        acb_sub_ui(t, s, 1, prec);
+        acb_get_mag(r, t);
+        mag_mul_ui(u, r, d);
+
+        if (mag_cmp_2exp_si(u, 3) < 0)
+        {
+            acb_one(t);
+            _acb_poly_zeta_cpx_series(z, t, a, 1, d, prec);
+
+            for (n = 0; n < d; n++)
+            {
+                mag_mul_ui(u, r, n + 1);
+                mag_one(v);
+                mag_mul_2exp_si(v, v, 3);
+                mag_sub_lower(v, v, u);
+                mag_div(v, u, v);
+                mag_mul_2exp_si(v, v, -3 * n);
+
+                arb_add_error_mag(acb_realref(z + n), v);
+                if (!is_real)
+                    arb_add_error_mag(acb_imagref(z + n), v);
+            }
+        }
+        else
+        {
+            _acb_vec_indeterminate(z, d);
+        }
+
+        mag_clear(r);
+        mag_clear(u);
+        mag_clear(v);
+        acb_clear(t);
+
+        return;
+    }
+
     is_real = const_is_real = 0;
 
     if (acb_is_real(s) && acb_is_real(a))