first implementation of Lerch transcendent

acb_dirichlet.h

@@ -48,6 +48,10 @@ void acb_dirichlet_zeta_jet(acb_t res, const acb_t s, int deflate, slong len, sl
 
 void acb_dirichlet_hurwitz(acb_t res, const acb_t s, const acb_t a, slong prec);
 
+void acb_dirichlet_lerch_phi_integral(acb_t res, const acb_t z, const acb_t s, const acb_t a, slong prec);
+void acb_dirichlet_lerch_phi_direct(acb_t res, const acb_t z, const acb_t s, const acb_t a, slong prec);
+void acb_dirichlet_lerch_phi(acb_t res, const acb_t z, const acb_t s, const acb_t a, slong prec);
+
 void acb_dirichlet_stieltjes(acb_t res, const fmpz_t n, const acb_t a, slong prec);
 
 typedef struct

acb_dirichlet/lerch_phi.c

@@ -0,0 +1,153 @@
+/*
+    Copyright (C) 2022 Fredrik Johansson
+
+    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_hypgeom.h"
+#include "acb_dirichlet.h"
+
+void
+acb_dirichlet_lerch_phi(acb_t res, const acb_t z, const acb_t s, const acb_t a, slong prec)
+{
+    if (!acb_is_finite(z) || !acb_is_finite(s) || !acb_is_finite(a))
+    {
+        acb_indeterminate(res);
+        return;
+    }
+
+    if (acb_contains_int(a) && !arb_is_positive(acb_realref(a)))
+    {
+        if (!(acb_is_int(s) && arb_is_nonpositive(acb_realref(s))))
+        {
+            acb_indeterminate(res);
+            return;
+        }
+    }
+
+    if (acb_is_zero(z))
+    {
+        acb_t t;
+        acb_init(t);
+        acb_neg(t, s);
+        acb_pow(res, a, t, prec);
+        acb_clear(t);
+        return;
+    }
+
+    if (acb_is_one(z))
+    {
+        arb_t one;
+        arb_init(one);
+        if (arb_gt(acb_realref(s), one))
+            acb_dirichlet_hurwitz(res, s, a, prec);
+        else
+            acb_indeterminate(res);
+        arb_clear(one);
+        return;
+    }
+
+    if (acb_equal_si(z, -1))
+    {
+        if (acb_is_one(a))
+        {
+            acb_dirichlet_eta(res, s, prec);
+        }
+        else if (acb_is_one(s))
+        {
+            /* (psi((a+1)/2) - psi(a/2))/2 */
+            acb_t t, u;
+            acb_init(t);
+            acb_init(u);
+            acb_mul_2exp_si(t, a, -1);
+            acb_digamma(t, t, prec);
+            acb_add_ui(u, a, 1, prec);
+            acb_mul_2exp_si(u, u, -1);
+            acb_digamma(u, u, prec);
+            acb_sub(res, u, t, prec);
+            acb_mul_2exp_si(res, res, -1);
+            acb_clear(t);
+            acb_clear(u);
+        }
+        else
+        {
+            /* 2^(-s) (zeta(s,a/2) - zeta(s,(a+1)/2)) */
+            acb_t t, u;
+            acb_init(t);
+            acb_init(u);
+            acb_mul_2exp_si(t, a, -1);
+            acb_hurwitz_zeta(t, s, t, prec);
+            acb_add_ui(u, a, 1, prec);
+            acb_mul_2exp_si(u, u, -1);
+            acb_hurwitz_zeta(u, s, u, prec);
+            acb_sub(t, t, u, prec);
+            acb_neg(u, s);
+            acb_set_ui(res, 2);
+            acb_pow(res, res, u, prec);
+            acb_mul(res, res, t, prec);
+            acb_clear(t);
+            acb_clear(u);
+        }
+        return;
+    }
+
+    if (acb_is_zero(s))
+    {
+        acb_sub_ui(res, z, 1, prec + 5);
+        acb_neg(res, res);
+        acb_inv(res, res, prec);
+        return;
+    }
+
+    if (acb_is_one(s))
+    {
+        acb_t t, u;
+        acb_init(t);
+        acb_init(u);
+        acb_one(t);
+        acb_add_ui(u, a, 1, prec + 5);
+        acb_hypgeom_2f1(t, t, a, u, z, ACB_HYPGEOM_2F1_BC, prec + 5);
+        acb_div(res, t, a, prec);
+        if (!acb_is_finite(res))
+            acb_indeterminate(res);
+        acb_clear(t);
+        acb_clear(u);
+        return;
+    }
+
+    if (acb_is_one(a) && !acb_contains_zero(z))
+    {
+        acb_t t;
+        acb_init(t);
+        acb_polylog(t, s, z, prec);
+        acb_div(res, t, z, prec);
+        acb_clear(t);
+        return;
+    }
+
+    {
+        mag_t zm, lim;
+        mag_init(zm);
+        mag_init(lim);
+
+        acb_get_mag(zm, z);
+        mag_set_d(lim, 0.75);
+
+        if (mag_cmp(zm, lim) <= 0)
+        {
+            acb_dirichlet_lerch_phi_direct(res, z, s, a, prec);
+        }
+        else
+        {
+            acb_dirichlet_lerch_phi_integral(res, z, s, a, prec);
+        }
+
+        mag_clear(zm);
+        mag_clear(lim);
+    }
+}

acb_dirichlet/lerch_phi_direct.c

@@ -0,0 +1,150 @@
+/*
+    Copyright (C) 2022 Fredrik Johansson
+
+    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_dirichlet_lerch_phi_direct(acb_t res, const acb_t z, const acb_t s, const acb_t a, slong prec)
+{
+    slong N, Nmax, wp, n;
+    int a_real;
+    acb_t negs, t, u, sum;
+    mag_t C, S, zmag, tail_bound, tm, tol;
+
+    if (!acb_is_finite(z) || !acb_is_finite(s) || !acb_is_finite(a))
+    {
+        acb_indeterminate(res);
+        return;
+    }
+
+    if (acb_contains_int(a) && !arb_is_positive(acb_realref(a)))
+    {
+        if (!(acb_is_int(s) && arb_is_nonpositive(acb_realref(s))))
+        {
+            acb_indeterminate(res);
+            return;
+        }
+    }
+
+    acb_init(negs);
+    acb_init(t);
+    acb_init(u);
+    acb_init(sum);
+
+    acb_neg(negs, s);
+
+    mag_init(C);
+    mag_init(S);
+    mag_init(zmag);
+    mag_init(tail_bound);
+    mag_init(tm);
+    mag_init(tol);
+
+    a_real = acb_is_real(a);
+    wp = prec + 10;
+
+    acb_get_mag(zmag, z);
+
+    /* first term: 1/a^s */
+    acb_pow(sum, a, negs, wp);
+
+    acb_get_mag(tol, sum);
+    mag_mul_2exp_si(tol, tol, -wp);
+
+    if (a_real)
+    {
+        /* Tail bound |z|^N / |(a+N)^s| * sum C^k, C = |z| * exp(max(0, -re(s)) / (a+N))  */
+        arb_nonnegative_part(acb_realref(t), acb_realref(negs));
+        arb_get_mag(S, acb_realref(t));
+    }
+    else
+    {
+        /* Tail bound |z|^N / |(a+N)^s| * sum C^k, C = |z| * exp(|s / (a+N)|)  */
+        acb_get_mag(S, s);
+    }
+
+    Nmax = 100 * prec + 0.1 * prec * n_sqrt(prec);
+    Nmax = FLINT_MAX(Nmax, 1);
+    Nmax = FLINT_MIN(Nmax, WORD_MAX / 2);
+
+    mag_inf(tail_bound);
+
+    for (N = 1; N <= Nmax; N = FLINT_MAX(N+4, N*1.1))
+    {
+        acb_add_ui(t, a, N, 53);
+
+        if (arb_is_positive(acb_realref(t)))
+        {
+            acb_get_mag_lower(C, t);
+            mag_div(C, S, C);
+            mag_exp(C, C);
+            mag_mul(C, C, zmag);
+            mag_geom_series(C, C, 0);
+
+            if (mag_is_finite(C))
+            {
+                mag_pow_ui(tail_bound, zmag, N);
+                mag_mul(tail_bound, tail_bound, C);
+                acb_pow(t, t, negs, 53);
+                acb_get_mag(C, t);
+                mag_mul(tail_bound, tail_bound, C);
+
+                if (mag_cmp(tail_bound, tol) <= 0)
+                    break;
+            }
+            else
+            {
+                mag_inf(tail_bound);
+            }
+        }
+    }
+
+    if (mag_is_finite(tail_bound))
+    {
+        acb_one(t);
+
+        for (n = 1; n < N; n++)
+        {
+            if (n % 8 == 0 && !acb_is_real(z))
+                acb_pow_ui(t, z, n, wp);
+            else
+                acb_mul(t, t, z, wp);
+
+            acb_add_ui(u, a, n, wp);
+            acb_pow(u, u, negs, wp);
+            acb_mul(u, t, u, wp);
+            acb_add(sum, sum, u, wp);
+        }
+
+        if (acb_is_real(z) && acb_is_real(s) && acb_is_real(a))
+            arb_add_error_mag(acb_realref(sum), tail_bound);
+        else
+            acb_add_error_mag(sum, tail_bound);
+
+        acb_set_round(res, sum, prec);
+    }
+    else
+    {
+        acb_indeterminate(res);
+    }
+
+    mag_clear(C);
+    mag_clear(S);
+    mag_clear(zmag);
+    mag_clear(tail_bound);
+    mag_clear(tm);
+    mag_clear(tol);
+
+    acb_clear(negs);
+    acb_clear(t);
+    acb_clear(u);
+    acb_clear(sum);
+}