riemann xi changes

acb_dirichlet.h

@@ -45,8 +45,6 @@ void acb_dirichlet_zeta(acb_t res, const acb_t s, slong prec);
 void acb_dirichlet_zeta_jet_rs(acb_ptr res, const acb_t s, slong len, slong prec);
 void acb_dirichlet_zeta_jet(acb_t res, const acb_t s, int deflate, slong len, slong prec);
 
-void acb_riemann_xi(acb_t res, const acb_t s, slong prec);
-
 void acb_dirichlet_hurwitz(acb_t res, const acb_t s, const acb_t a, slong prec);
 
 typedef struct
@@ -75,6 +73,8 @@ void _acb_dirichlet_euler_product_real_ui(arb_t res, ulong s,
 
 void acb_dirichlet_eta(acb_t res, const acb_t s, slong prec);
 
+void acb_dirichlet_xi(acb_t res, const acb_t s, slong prec);
+
 void acb_dirichlet_pairing(acb_t res, const dirichlet_group_t G, ulong m, ulong n, slong prec);
 void acb_dirichlet_pairing_char(acb_t res, const dirichlet_group_t G, const dirichlet_char_t a, const dirichlet_char_t b, slong prec);
 

acb_dirichlet/riemann_xi.c

@@ -1,51 +0,0 @@
-/*
-    Copyright (C) 2018 D.H.J Polymath
-
-    This file is part of Arb.
-
-    Arb is free software: you can redistribute it and/or modify it under
-    the terms of the GNU Lesser General Public License (LGPL) as published
-    by the Free Software Foundation; either version 2.1 of the License, or
-    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
-*/
-
-#include "acb_dirichlet.h"
-
-void acb_riemann_xi(acb_t res, const acb_t s, slong prec)
-{
-    acb_t pi, z1, z2, z3, z4;
-
-    acb_init(pi);
-    acb_const_pi(pi, prec);
-
-    /* s*(s-1)/2 */
-    acb_init(z1);
-    acb_sub_ui(z1, s, 1, prec);
-    acb_mul(z1, z1, s, prec);
-    acb_mul_2exp_si(z1, z1, -1);
-
-    /* pi^(-s/2) */
-    acb_init(z2);
-    acb_mul_2exp_si(z2, s, -1);
-    acb_neg(z2, z2);
-    acb_pow(z2, pi, z2, prec);
-
-    /* gamma(s/2) */
-    acb_init(z3);
-    acb_mul_2exp_si(z3, s, -1);
-    acb_gamma(z3, z3, prec);
-
-    /* zeta(s) */
-    acb_init(z4);
-    acb_zeta(z4, s, prec);
-
-    acb_mul(res, z1, z2, prec);
-    acb_mul(res, res, z3, prec);
-    acb_mul(res, res, z4, prec);
-
-    acb_clear(pi);
-    acb_clear(z1);
-    acb_clear(z2);
-    acb_clear(z3);
-    acb_clear(z4);
-}

acb_dirichlet/test/t-xi.c

@@ -16,7 +16,7 @@ int main()
     slong iter;
     flint_rand_t state;
 
-    flint_printf("riemann_xi....");
+    flint_printf("xi....");
     fflush(stdout);
 
     flint_randinit(state);
@@ -36,8 +36,8 @@ int main()
         arb_randtest_precise(acb_realref(a), state, 1 + n_randint(state, 1000), 3);
         arb_randtest_precise(acb_imagref(a), state, 1 + n_randint(state, 1000), 3);
 
-        acb_riemann_xi(b, a, prec1);
-        acb_riemann_xi(c, a, prec2);
+        acb_dirichlet_xi(b, a, prec1);
+        acb_dirichlet_xi(c, a, prec2);
 
         if (!acb_overlaps(b, c))
         {
@@ -48,10 +48,22 @@ int main()
             flint_abort();
         }
 
-        /* check riemann_xi(s) = riemann_xi(1-s) */
+        acb_set(c, a);
+        acb_dirichlet_xi(c, c, prec1);
+
+        if (!acb_equal(b, c))
+        {
+            flint_printf("FAIL: aliasing\n\n");
+            flint_printf("a = "); acb_printd(a, 30); flint_printf("\n\n");
+            flint_printf("b = "); acb_printd(b, 30); flint_printf("\n\n");
+            flint_printf("c = "); acb_printd(c, 30); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        /* check xi(s) = xi(1-s) */
         acb_sub_ui(c, a, 1, prec1);
         acb_neg(c, c);
-        acb_riemann_xi(c, c, prec1);
+        acb_dirichlet_xi(c, c, prec1);
 
         if (!acb_overlaps(b, c))
         {

acb_dirichlet/xi.c

@@ -0,0 +1,68 @@
+/*
+    Copyright (C) 2018 D.H.J Polymath
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_dirichlet.h"
+
+static void
+_acb_dirichlet_xi(acb_t res, const acb_t s, slong prec)
+{
+    acb_t pi, z1, z2, z3;
+
+    acb_init(pi);
+    acb_init(z1);
+    acb_init(z2);
+    acb_init(z3);
+
+    /* (s-1) * pi^(-s/2) * gamma(1 + s/2) * zeta(s) */
+    acb_const_pi(pi, prec);
+    acb_mul_2exp_si(z1, s, -1);
+    acb_neg(z1, z1);
+    acb_pow(z1, pi, z1, prec);
+    acb_mul_2exp_si(z2, s, -1);
+    acb_add_ui(z2, z2, 1, prec);
+    acb_gamma(z2, z2, prec);
+    acb_zeta(z3, s, prec);
+    acb_sub_ui(res, s, 1, prec);
+    acb_mul(res, res, z1, prec);
+    acb_mul(res, res, z2, prec);
+    acb_mul(res, res, z3, prec);
+
+    acb_clear(pi);
+    acb_clear(z1);
+    acb_clear(z2);
+    acb_clear(z3);
+}
+
+/* todo: interval containing 1 */
+void acb_dirichlet_xi(acb_t res, const acb_t s, slong prec)
+{
+    if (!acb_is_finite(s))
+    {
+        acb_indeterminate(res);
+    }
+    else if (acb_is_one(s))
+    {
+        acb_one(res);
+        acb_mul_2exp_si(res, res, -1);
+    }
+    else if (arf_sgn(arb_midref(acb_realref(s))) < 0)
+    {
+        acb_t t;
+        acb_sub_ui(t, s, 1, prec);
+        acb_neg(t, t);
+        _acb_dirichlet_xi(res, t, prec);
+        acb_clear(t);
+    }
+    else
+    {
+        _acb_dirichlet_xi(res, s, prec);
+    }
+}

doc/source/acb_dirichlet.rst

@@ -123,6 +123,12 @@ Riemann zeta function
     Note that the alternating character `\{1,-1\}` is not itself
     a Dirichlet character.
 
+.. function:: void acb_dirichlet_xi(acb_t res, const acb_t s, slong prec)
+
+    Sets *res* to the Riemann xi function
+    `\xi(s) = \frac{1}{2} s (s-1) \pi^{-s/2} \Gamma(\frac{1}{2} s) \zeta(s)`.
+    The functional equation for xi is `\xi(1-s) = \xi(s)`.
+
 Riemann-Siegel formula
 -------------------------------------------------------------------------------