.. _acb-dirichlet:

**acb_dirichlet.h** -- Dirichlet L-functions, Riemann zeta and related functions
===================================================================================

This module allows working with values of Dirichlet characters,
Dirichlet L-functions, and related functions.
A Dirichlet L-function is the analytic continuation of an L-series

.. math ::

    L(s,\chi) = \sum_{k=1}^\infty \frac{\chi(k)}{k^s}

where `\chi(k)` is a Dirichlet character. The trivial character
`\chi(k) = 1` gives the Riemann zeta function.
Working with Dirichlet characters is documented in :ref:`dirichlet`.

The code in other modules for computing the Riemann zeta function,
Hurwitz zeta function and polylogarithm will possibly be migrated to this
module in the future.

Roots of unity
-------------------------------------------------------------------------------

.. type:: acb_dirichlet_roots_struct

.. type:: acb_dirichlet_roots_t

.. function:: void acb_dirichlet_roots_init(acb_dirichlet_roots_t roots, ulong n, slong num, slong prec)

    Initializes *roots* with precomputed data for fast evaluation of roots of
    unity `e^{2\pi i k/n}` of a fixed order *n*. The precomputation is
    optimized for *num* evaluations.

    For very small *num*, only the single root `e^{2\pi i/n}` will be
    precomputed, which can then be raised to a power. For small *prec*
    and large *n*, this method might even skip precomputing this single root
    if it estimates that evaluating roots of unity from scratch will be faster
    than powering.

    If *num* is large enough, the whole set of roots in the first quadrant
    will be precomputed at once. However, this is automatically avoided for
    large *n* if too much memory would be used. For intermediate *num*,
    baby-step giant-step tables are computed.

.. function:: void acb_dirichlet_roots_clear(acb_dirichlet_roots_t roots)

    Clears the structure.

.. function:: void acb_dirichlet_root(acb_t res, const acb_dirichlet_roots_t roots, ulong k, slong prec)

    Computes `e^{2\pi i k/n}`.

Truncated L-series and power sums
-------------------------------------------------------------------------------

.. function:: void acb_dirichlet_powsum_term(acb_ptr res, arb_t log_prev, ulong * prev, const acb_t s, ulong k, int integer, int critical_line, slong len, slong prec)

    Sets *res* to `k^{-(s+x)}` as a power series in *x* truncated to length *len*.
    The flags *integer* and *critical_line* respectively specify optimizing
    for *s* being an integer or having real part 1/2.

    On input *log_prev* should contain the natural logarithm of the integer
    at *prev*. If *prev* is close to *k*, this can be used to speed up
    computations. If `\log(k)` is computed internally by this function, then
    *log_prev* is overwritten by this value, and the integer at *prev* is
    overwritten by *k*, allowing *log_prev* to be recycled for the next
    term when evaluating a power sum.

.. function:: void acb_dirichlet_powsum_sieved(acb_ptr res, const acb_t s, ulong n, slong len, slong prec)

    Sets *res* to `\sum_{k=1}^n k^{-(s+x)}`
    as a power series in *x* truncated to length *len*.
    This function stores a table of powers that have already been calculated,
    computing `(ij)^r` as `i^r j^r` whenever `k = ij` is
    composite. As a further optimization, it groups all even `k` and
    evaluates the sum as a polynomial in `2^{-(s+x)}`.
    This scheme requires about `n / \log n` powers, `n / 2` multiplications,
    and temporary storage of `n / 6` power series. Due to the extra
    power series multiplications, it is only faster than the naive
    algorithm when *len* is small.

.. function:: void acb_dirichlet_powsum_smooth(acb_ptr res, const acb_t s, ulong n, slong len, slong prec)

    Sets *res* to `\sum_{k=1}^n k^{-(s+x)}`
    as a power series in *x* truncated to length *len*.
    This function performs partial sieving by adding multiples of 5-smooth *k*
    into separate buckets. Asymptotically, this requires computing 4/15
    of the powers, which is slower than *sieved*, but only requires
    logarithmic extra space. It is also faster for large *len*, since most
    power series multiplications are traded for additions.
    A slightly bigger gain for larger *n* could be achieved by using more
    small prime factors, at the expense of space.

Riemann zeta function
-------------------------------------------------------------------------------

.. function:: void acb_dirichlet_zeta(acb_t res, const acb_t s, slong prec)

    Computes `\zeta(s)` using an automatic choice of algorithm.

.. function:: void acb_dirichlet_zeta_jet(acb_t res, const acb_t s, int deflate, slong len, slong prec)

    Computes the first *len* terms of the Taylor series of the Riemann zeta
    function at *s*. If *deflate* is nonzero, computes the deflated
    function `\zeta(s) - 1/(s-1)` instead.

.. function:: void acb_dirichlet_zeta_bound(mag_t res, const acb_t s)

    Computes an upper bound for `|\zeta(s)|` quickly. On the critical strip (and
    slightly outside of it), formula (43.3) in [Rad1973]_ is used.
    To the right, evaluating at the real part of *s* gives a trivial bound.
    To the left, the functional equation is used.

.. function:: void acb_dirichlet_zeta_deriv_bound(mag_t der1, mag_t der2, const acb_t s)

    Sets *der1* to a bound for `|\zeta'(s)|` and *der2* to a bound for
    `|\zeta''(s)|`. These bounds are mainly intended for use in the critical
    strip and will not be tight.

.. function:: void acb_dirichlet_eta(acb_t res, const acb_t s, slong prec)

    Sets *res* to the Dirichlet eta function
    `\eta(s) = \sum_{k=1}^{\infty} (-1)^{k+1} / k^s = (1-2^{1-s}) \zeta(s)`,
    also known as the alternating 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
-------------------------------------------------------------------------------

The Riemann-Siegel (RS) formula is implemented closely following
J. Arias de Reyna [Ari2011]_.
For `s = \sigma + it` with `t > 0`, the expansion takes the form

.. math ::

    \zeta(s) = \mathcal{R}(s) + X(s) \overline{\mathcal{R}}(1-s), \quad X(s) = \pi^{s-1/2} \frac{\Gamma((1-s)/2)}{\Gamma(s/2)}

where

.. math ::

    \mathcal{R}(s) = \sum_{k=1}^N \frac{1}{k^s} + (-1)^{N-1} U a^{-\sigma} \left[ \sum_{k=0}^K \frac{C_k(p)}{a^k} + RS_K \right]

.. math ::

    U = \exp\left(-i\left[ \frac{t}{2} \log\left(\frac{t}{2\pi}\right)-\frac{t}{2}-\frac{\pi}{8} \right]\right), \quad
    a = \sqrt{\frac{t}{2\pi}}, \quad N = \lfloor a \rfloor, \quad p = 1-2(a-N).

The coefficients `C_k(p)` in the asymptotic part of the expansion
are expressed in terms of certain auxiliary coefficients `d_j^{(k)}`
and `F^{(j)}(p)`.
Because of artificial discontinuities, *s* should be exact inside
the evaluation.

.. function:: void acb_dirichlet_zeta_rs_f_coeffs(acb_ptr f, const arb_t p, slong n, slong prec)

    Computes the coefficients `F^{(j)}(p)` for `0 \le j < n`.
    Uses power series division. This method breaks down when `p = \pm 1/2`
    (which is not problem if *s* is an exact floating-point number).

.. function:: void acb_dirichlet_zeta_rs_d_coeffs(arb_ptr d, const arb_t sigma, slong k, slong prec)

    Computes the coefficients `d_j^{(k)}` for `0 \le j \le \lfloor 3k/2 \rfloor + 1`.
    On input, the array *d* must contain the coefficients for `d_j^{(k-1)}`
    unless `k = 0`, and these coefficients will be updated in-place.

.. function:: void acb_dirichlet_zeta_rs_bound(mag_t err, const acb_t s, slong K)

    Bounds the error term `RS_K` following Theorem 4.2 in Arias de Reyna.

.. function:: void acb_dirichlet_zeta_rs_r(acb_t res, const acb_t s, slong K, slong prec)

    Computes `\mathcal{R}(s)` in the upper half plane. Uses precisely *K*
    asymptotic terms in the RS formula if this input parameter is positive;
    otherwise chooses the number of terms automatically based on *s* and the
    precision.

.. function:: void acb_dirichlet_zeta_rs(acb_t res, const acb_t s, slong K, slong prec)

    Computes `\zeta(s)` using the Riemann-Siegel formula. Uses precisely
    *K* asymptotic terms in the RS formula if this input parameter is positive;
    otherwise chooses the number of terms automatically based on *s* and the
    precision.

.. function:: void acb_dirichlet_zeta_jet_rs(acb_t res, const acb_t s, slong len, slong prec)

    Computes the first *len* terms of the Taylor series of the Riemann zeta
    function at *s* using the Riemann Siegel formula. This function currently
    only supports *len* = 1 or *len* = 2. A finite difference is used
    to compute the first derivative.

Hurwitz zeta function
-------------------------------------------------------------------------------

.. function:: void acb_dirichlet_hurwitz(acb_t res, const acb_t s, const acb_t a, slong prec)

    Computes the Hurwitz zeta function `\zeta(s, a)`.
    This function automatically delegates to the code for the Riemann zeta function
    when `a = 1`. Some other special cases may also be handled by direct
    formulas. In general, Euler-Maclaurin summation is used.

Hurwitz zeta function precomputation
-------------------------------------------------------------------------------

.. type:: acb_dirichlet_hurwitz_precomp_struct

.. type:: acb_dirichlet_hurwitz_precomp_t

.. function:: void acb_dirichlet_hurwitz_precomp_init(acb_dirichlet_hurwitz_precomp_t pre, const acb_t s, int deflate, ulong A, ulong K, ulong N, slong prec)

    Precomputes a grid of Taylor polynomials for fast evaluation of
    `\zeta(s,a)` on `a \in (0,1]` with fixed *s*.
    *A* is the initial shift to apply to *a*, *K* is the number of Taylor terms,
    *N* is the number of grid points.  The precomputation requires *NK*
    evaluations of the Hurwitz zeta function, and each subsequent evaluation
    requires *2K* simple arithmetic operations (polynomial evaluation) plus
    *A* powers. As *K* grows, the error is at most `O(1/(2AN)^K)`.

    This function can be called with *A* set to zero, in which case
    no Taylor series precomputation is performed. This means that evaluation
    will be identical to calling :func:`acb_dirichlet_hurwitz` directly.

    Otherwise, we require that *A*, *K* and *N* are all positive. For a finite
    error bound, we require `K+\operatorname{re}(s) > 1`.
    To avoid an initial "bump" that steals precision
    and slows convergence, *AN* should be at least roughly as large as `|s|`,
    e.g. it is a good idea to have at least `AN > 0.5 |s|`.

    If *deflate* is set, the deflated Hurwitz zeta function is used,
    removing the pole at `s = 1`.

.. function:: void acb_dirichlet_hurwitz_precomp_init_num(acb_dirichlet_hurwitz_precomp_t pre, const acb_t s, int deflate, double num_eval, slong prec)

    Initializes *pre*, choosing the parameters *A*, *K*, and *N*
    automatically to minimize the cost of *num_eval* evaluations of the
    Hurwitz zeta function at argument *s* to precision *prec*.

.. function:: void acb_dirichlet_hurwitz_precomp_clear(acb_dirichlet_hurwitz_precomp_t pre)

    Clears the precomputed data.

.. function:: void acb_dirichlet_hurwitz_precomp_choose_param(ulong * A, ulong * K, ulong * N, const acb_t s, double num_eval, slong prec)

    Chooses precomputation parameters *A*, *K* and *N* to minimize
    the cost of *num_eval* evaluations of the Hurwitz zeta function
    at argument *s* to precision *prec*.
    If it is estimated that evaluating each Hurwitz zeta function from
    scratch would be better than performing a precomputation, *A*, *K* and *N*
    are all set to 0.

.. function:: void acb_dirichlet_hurwitz_precomp_bound(mag_t res, const acb_t s, ulong A, ulong K, ulong N)

    Computes an upper bound for the truncation error (not accounting for
    roundoff error) when evaluating `\zeta(s,a)` with precomputation parameters
    *A*, *K*, *N*, assuming that `0 < a \le 1`.
    For details, see :ref:`algorithms_hurwitz`.

.. function:: void acb_dirichlet_hurwitz_precomp_eval(acb_t res, const acb_dirichlet_hurwitz_precomp_t pre, ulong p, ulong q, slong prec)

    Evaluates `\zeta(s,p/q)` using precomputed data, assuming that `0 < p/q \le 1`.

Stieltjes constants
-------------------------------------------------------------------------------

.. function:: void acb_dirichlet_stieltjes(acb_t res, const fmpz_t n, const acb_t a, slong prec)

    Given a nonnegative integer *n*, sets *res* to the generalized Stieltjes constant
    `\gamma_n(a)` which is the coefficient in the Laurent series of the
    Hurwitz zeta function at the pole

    .. math ::

        \zeta(s,a) = \frac{1}{s-1} + \sum_{n=0}^\infty \frac{(-1)^n}{n!} \gamma_n(a) (s-1)^n.

    With `a = 1`, this gives the ordinary Stieltjes constants for the
    Riemann zeta function.

    This function uses an integral representation to permit fast computation
    for extremely large *n* [JB2018]_. If *n* is moderate and the precision
    is high enough, it falls back to evaluating the Hurwitz zeta function
    of a power series and reading off the last coefficient.

    Note that for computing a range of values
    `\gamma_0(a), \ldots, \gamma_n(a)`, it is
    generally more efficient to evaluate the Hurwitz zeta function series
    expansion once at `s = 1` than to call this function repeatedly,
    unless *n* is extremely large (at least several hundred).

Dirichlet character evaluation
-------------------------------------------------------------------------------

.. function:: void acb_dirichlet_chi(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, ulong n, slong prec)

    Sets *res* to `\chi(n)`, the value of the Dirichlet character *chi*
    at the integer *n*.

.. function:: void acb_dirichlet_chi_vec(acb_ptr v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv, slong prec)

    Compute the *nv* first Dirichlet values.

.. function:: void acb_dirichlet_pairing(acb_t res, const dirichlet_group_t G, ulong m, ulong n, slong prec)

.. function:: 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)

    Sets *res* to the value of the Dirichlet pairing `\chi(m,n)` at numbers `m` and `n`.
    The second form takes two characters as input.

Dirichlet character Gauss, Jacobi and theta sums
-------------------------------------------------------------------------------

.. function:: void acb_dirichlet_gauss_sum_naive(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)

.. function:: void acb_dirichlet_gauss_sum_factor(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)

.. function:: void acb_dirichlet_gauss_sum_order2(acb_t res, const dirichlet_char_t chi, slong prec)

.. function:: void acb_dirichlet_gauss_sum_theta(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)

.. function:: void acb_dirichlet_gauss_sum(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)

.. function:: void acb_dirichlet_gauss_sum_ui(acb_t res, const dirichlet_group_t G, ulong a, slong prec)

   Sets *res* to the Gauss sum

   .. math::

      G_q(a) = \sum_{x \bmod q} \chi_q(a, x) e^{\frac{2i\pi x}q}

   - the *naive* version computes the sum as defined.

   - the *factor* version writes it as a product of local Gauss sums by chinese
     remainder theorem.

   - the *order2* version assumes *chi* is real and primitive and returns
     `i^p\sqrt q` where `p` is the parity of `\chi`.

   - the *theta* version assumes that *chi* is primitive to obtain the Gauss
     sum by functional equation of the theta series at `t=1`. An abort will be
     raised if the theta series vanishes at `t=1`. Only 4 exceptional
     characters of conductor 300 and 600 are known to have this particularity,
     and none with primepower modulus.

   - the default version automatically combines the above methods.

   - the *ui* version only takes the Conrey number *a* as parameter.

.. function:: void acb_dirichlet_jacobi_sum_naive(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi1, const dirichlet_char_t chi2, slong prec)

.. function:: void acb_dirichlet_jacobi_sum_factor(acb_t res,  const dirichlet_group_t G, const dirichlet_char_t chi1, const dirichlet_char_t chi2, slong prec)

.. function:: void acb_dirichlet_jacobi_sum_gauss(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi1, const dirichlet_char_t chi2, slong prec)

.. function:: void acb_dirichlet_jacobi_sum(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi1,  const dirichlet_char_t chi2, slong prec)

.. function:: void acb_dirichlet_jacobi_sum_ui(acb_t res, const dirichlet_group_t G, ulong a, ulong b, slong prec)

   Computes the Jacobi sum

   .. math::

      J_q(a,b) = \sum_{x \bmod q} \chi_q(a, x)\chi_q(b, 1-x)

   - the *naive* version computes the sum as defined.

   - the *factor* version writes it as a product of local Jacobi sums

   - the *gauss* version assumes `ab` is primitive and uses the formula
     `J_q(a,b)G_q(ab) = G_q(a)G_q(b)`

   - the default version automatically combines the above methods.

   - the *ui* version only takes the Conrey numbers *a* and *b* as parameters.

.. function:: void acb_dirichlet_chi_theta_arb(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, const arb_t t, slong prec)

.. function:: void acb_dirichlet_ui_theta_arb(acb_t res, const dirichlet_group_t G, ulong a, const arb_t t, slong prec)

    Compute the theta series `\Theta_q(a,t)` for real argument `t>0`.
    Beware that if `t<1` the functional equation

    .. math::

        t \theta(a,t) = \epsilon(\chi) \theta\left(\frac1a, \frac1t\right)

    should be used, which is not done automatically (to avoid recomputing the
    Gauss sum).

    We call *theta series* of a Dirichlet character the quadratic series

    .. math::

       \Theta_q(a) = \sum_{n\geq 0} \chi_q(a, n) n^p x^{n^2}

    where `p` is the parity of the character `\chi_q(a,\cdot)`.

    For `\Re(t)>0` we write `x(t)=\exp(-\frac{\pi}{N}t^2)` and define

    .. math::

       \Theta_q(a,t) = \sum_{n\geq 0} \chi_q(a, n) x(t)^{n^2}.

.. function:: ulong acb_dirichlet_theta_length(ulong q, const arb_t t, slong prec)

   Compute the number of terms to be summed in the theta series of argument *t*
   so that the tail is less than `2^{-\mathrm{prec}}`.

.. function:: void acb_dirichlet_qseries_powers_naive(acb_t res, const arb_t x, int p, const ulong * a, const acb_dirichlet_powers_t z, slong len, slong prec)

.. function:: void acb_dirichlet_qseries_powers_smallorder(acb_t res, const arb_t x, int p, const ulong * a, const acb_dirichlet_powers_t z, slong len, slong prec)

   Compute the series `\sum n^p z^{a_n} x^{n^2}` for exponent list *a*,
   precomputed powers *z* and parity *p* (being 0 or 1).

   The *naive* version sums the series as defined, while the *smallorder*
   variant evaluates the series on the quotient ring by a cyclotomic polynomial
   before evaluating at the root of unity, ignoring its argument *z*.

Discrete Fourier transforms
-------------------------------------------------------------------------------

If `f` is a function `\mathbb Z/q\mathbb Z\to \mathbb C`,
its discrete Fourier transform is the function
defined on Dirichlet characters mod `q` by

.. math::

   \hat f(\chi) = \sum_{x\mod q}\overline{\chi(x)}f(x)

See the :ref:`acb-dft` module.

Here we take advantage of the Conrey isomorphism `G \to \hat G`
to consider the Fourier transform on Conrey labels as

.. math::

   g(a) = \sum_{b\bmod q}\overline{\chi_q(a,b)}f(b)


.. function:: void acb_dirichlet_dft_conrey(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec)

   Compute the DFT of *v* using Conrey indices.
   This function assumes *v* and *w* are vectors
   of size *G->phi_q*, whose values correspond to a lexicographic ordering
   of Conrey logs (as obtained using :func:`dirichlet_char_next` or
   by :func:`dirichlet_char_index`).

   For example, if `q=15`, the Conrey elements are stored in following
   order

   =======  =============  =====================
    index    log = [e,f]     number = 7^e 11^f
   =======  =============  =====================
      0       [0, 0]        1
      1       [0, 1]        7
      2       [0, 2]        4
      3       [0, 3]        13
      4       [0, 4]        1
      5       [1, 0]        11
      6       [1, 1]        2
      7       [1, 2]        14
      8       [1, 3]        8
      9       [1, 4]        11
   =======  =============  =====================

.. function:: void acb_dirichlet_dft(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec)

   Compute the DFT of *v* using Conrey numbers.
   This function assumes *v* and *w* are vectors of size *G->q*.
   All values at index not coprime to *G->q* are ignored.

Dirichlet L-functions
-------------------------------------------------------------------------------

.. function:: void acb_dirichlet_root_number_theta(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)

.. function:: void acb_dirichlet_root_number(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)

   Sets *res* to the root number `\epsilon(\chi)` for a primitive character *chi*,
   which appears in the functional equation (where `p` is the parity of `\chi`):

   .. math::

      \left(\frac{q}{\pi}\right)^{\frac{s+p}2}\Gamma\left(\frac{s+p}2\right) L(s, \chi) = \epsilon(\chi) \left(\frac{q}{\pi}\right)^{\frac{1-s+p}2}\Gamma\left(\frac{1-s+p}2\right) L(1 - s, \overline\chi)

   - The *theta* variant uses the evaluation at `t=1` of the Theta series.

   - The default version computes it via the gauss sum.

.. function:: void acb_dirichlet_l_hurwitz(acb_t res, const acb_t s, const acb_dirichlet_hurwitz_precomp_t precomp, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)

    Computes `L(s,\chi)` using decomposition in terms of the Hurwitz zeta function

    .. math::

        L(s,\chi) = q^{-s}\sum_{k=1}^q \chi(k) \,\zeta\!\left(s,\frac kq\right).

    If `s = 1` and `\chi` is non-principal, the deflated Hurwitz zeta function
    is used to avoid poles.

    If *precomp* is *NULL*, each Hurwitz zeta function value is computed
    directly. If a pre-initialized *precomp* object is provided, this will be
    used instead to evaluate the Hurwitz zeta function.

.. function:: void acb_dirichlet_l_euler_product(acb_t res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)

.. function:: void _acb_dirichlet_euler_product_real_ui(arb_t res, ulong s, const signed char * chi, int mod, int reciprocal, slong prec)

    Computes `L(s,\chi)` directly using the Euler product. This is
    efficient if *s* has large positive real part. As implemented, this
    function only gives a finite result if `\operatorname{re}(s) \ge 2`.

    An error bound is computed via :func:`mag_hurwitz_zeta_uiui`.
    If *s* is complex, replace it with its real part. Since

    .. math ::

        \frac{1}{L(s,\chi)} = \prod_{p} \left(1 - \frac{\chi(p)}{p^s}\right)
                = \sum_{k=1}^{\infty} \frac{\mu(k)\chi(k)}{k^s}

    and the truncated product gives all smooth-index terms in the series, we have

    .. math ::

        \left|\prod_{p < N} \left(1 - \frac{\chi(p)}{p^s}\right) - \frac{1}{L(s,\chi)}\right|
        \le \sum_{k=N}^{\infty} \frac{1}{k^s} = \zeta(s,N).

    The underscore version specialized for integer *s* assumes that `\chi` is
    a real Dirichlet character given by the explicit list *chi* of character
    values at 0, 1, ..., *mod* - 1. If *reciprocal* is set, it computes
    `1 / L(s,\chi)` (this is faster if the reciprocal can be used directly).

.. function:: void acb_dirichlet_l(acb_t res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)

    Computes `L(s,\chi)` using a default choice of algorithm.

.. function:: void acb_dirichlet_l_fmpq(acb_t res, const fmpq_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
              void acb_dirichlet_l_fmpq_afe(acb_t res, const fmpq_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)

    Computes `L(s,\chi)` where *s* is a rational number.
    The *afe* version uses the approximate functional equation;
    the default version chooses an algorithm automatically.

.. function:: void acb_dirichlet_l_vec_hurwitz(acb_ptr res, const acb_t s, const acb_dirichlet_hurwitz_precomp_t precomp, const dirichlet_group_t G, slong prec)

    Compute all values `L(s,\chi)` for `\chi` mod `q`, using the
    Hurwitz zeta function and a discrete Fourier transform.
    The output *res* is assumed to have length *G->phi_q* and values
    are stored by lexicographically ordered
    Conrey logs. See :func:`acb_dirichlet_dft_conrey`.

    If *precomp* is *NULL*, each Hurwitz zeta function value is computed
    directly. If a pre-initialized *precomp* object is provided, this will be
    used instead to evaluate the Hurwitz zeta function.

.. function:: void acb_dirichlet_l_jet(acb_ptr res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, int deflate, slong len, slong prec)

    Computes the Taylor expansion of `L(s,\chi)` to length *len*,
    i.e. `L(s), L'(s), \ldots, L^{(len-1)}(s) / (len-1)!`.
    If *deflate* is set, computes the expansion of

    .. math ::

        L(s,\chi) - \frac{\sum_{k=1}^q \chi(k)}{(s-1)q}

    instead. If *chi* is a principal character, then this has the effect of
    subtracting the pole with residue `\sum_{k=1}^q \chi(k) = \phi(q) / q`
    that is located at `s = 1`. In particular, when evaluated at `s = 1`, this
    gives the regular part of the Laurent expansion.
    When *chi* is non-principal, *deflate* has no effect.

.. function:: void _acb_dirichlet_l_series(acb_ptr res, acb_srcptr s, slong slen, const dirichlet_group_t G, const dirichlet_char_t chi, int deflate, slong len, slong prec)

.. function:: void acb_dirichlet_l_series(acb_poly_t res, const acb_poly_t s, const dirichlet_group_t G, const dirichlet_char_t chi, int deflate, slong len, slong prec)

    Sets *res* to the power series `L(s,\chi)` where *s* is a given power series, truncating the result to length *len*.
    See :func:`acb_dirichlet_l_jet` for the meaning of the *deflate* flag.

Hardy Z-functions
-------------------------------------------------------------------------------

For convenience, setting both *G* and *chi* to *NULL* in the following
methods selects the Riemann zeta function.

Currently, these methods require *chi* to be a primitive character.

.. function:: void acb_dirichlet_hardy_theta(acb_ptr res, const acb_t t, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)

    Computes the phase function used to construct the Z-function.
    We have

    .. math ::

        \theta(t) = -\frac{t}{2} \log(\pi/q) - \frac{i \log(\epsilon)}{2}
            + \frac{\log \Gamma((s+\delta)/2) - \log \Gamma((1-s+\delta)/2)}{2i}

    where `s = 1/2+it`, `\delta` is the parity of *chi*, and `\epsilon`
    is the root number as computed by :func:`acb_dirichlet_root_number`.
    The first *len* terms in the Taylor expansion are written to the output.

.. function:: void acb_dirichlet_hardy_z(acb_t res, const acb_t t, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)

    Computes the Hardy Z-function, also known as the Riemann-Siegel Z-function
    `Z(t) = e^{i \theta(t)} L(1/2+it)`, which is real-valued for real *t*.
    The first *len* terms in the Taylor expansion are written to the output.

.. function:: void _acb_dirichlet_hardy_theta_series(acb_ptr res, acb_srcptr t, slong tlen, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)

.. function:: void acb_dirichlet_hardy_theta_series(acb_poly_t res, const acb_poly_t t, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)

    Sets *res* to the power series `\theta(t)` where *t* is a given power series, truncating the result to length *len*.

.. function:: void _acb_dirichlet_hardy_z_series(acb_ptr res, acb_srcptr t, slong tlen, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)

.. function:: void acb_dirichlet_hardy_z_series(acb_poly_t res, const acb_poly_t t, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)

    Sets *res* to the power series `Z(t)` where *t* is a given power series, truncating the result to length *len*.

Gram points
-------------------------------------------------------------------------------

.. function:: void acb_dirichlet_gram_point(arb_t res, const fmpz_t n, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)

    Sets *res* to the *n*-th Gram point `g_n`, defined as the unique solution
    in `[7, \infty)` of `\theta(g_n) = \pi n`. Currently only the Gram points
    corresponding to the Riemann zeta function are supported and *G* and *chi*
    must both be set to *NULL*. Requires `n \ge -1`.

Riemann zeta function zeros
-------------------------------------------------------------------------------

The following functions for counting and isolating zeros of the Riemann zeta
function use the ideas from the implementation of Turing's method in
mpmath [Joh2018b]_ by Juan Arias de Reyna, described in [Ari2012]_.

.. function:: ulong acb_dirichlet_turing_method_bound(const fmpz_t p)

    Computes an upper bound *B* for the minimum number of consecutive good
    Gram blocks sufficient to count nontrivial zeros of the Riemann zeta
    function using Turing's method [Tur1953]_ as updated by [Leh1970]_,
    [Bre1979]_, and [Tru2011]_.

    Let `N(T)` denote the number of zeros (counted according to their
    multiplicities) of `\zeta(s)` in the region `0 < \operatorname{Im}(s) \le T`.
    If at least *B* consecutive Gram blocks with union `[g_n, g_p)`
    satisfy Rosser's rule, then `N(g_n) \le n + 1` and `N(g_p) \ge p + 1`.

.. function:: int _acb_dirichlet_definite_hardy_z(arb_t res, const arf_t t, slong * pprec)

    Sets *res* to the Hardy Z-function `Z(t)`.
    The initial precision (* *pprec*) is increased as necessary
    to determine the sign of `Z(t)`. The sign is returned.

.. function:: void _acb_dirichlet_isolate_gram_hardy_z_zero(arf_t a, arf_t b, const fmpz_t n)

    Uses Gram's law to compute an interval `(a, b)` that
    contains the *n*-th zero of the Hardy Z-function and no other zero.
    Requires `1 \le n \le 126`.

.. function:: void _acb_dirichlet_isolate_rosser_hardy_z_zero(arf_t a, arf_t b, const fmpz_t n)

    Uses Rosser's rule to compute an interval `(a, b)` that
    contains the *n*-th zero of the Hardy Z-function and no other zero.
    Requires `1 \le n \le 13999526`.

.. function:: void _acb_dirichlet_isolate_turing_hardy_z_zero(arf_t a, arf_t b, const fmpz_t n)

    Computes an interval `(a, b)` that contains the *n*-th zero of the
    Hardy Z-function and no other zero, following Turing's method.
    Requires `n \ge 2`.

.. function:: void acb_dirichlet_isolate_hardy_z_zero(arf_t a, arf_t b, const fmpz_t n)

    Computes an interval `(a, b)` that contains the *n*-th zero of the
    Hardy Z-function and contains no other zero, using the most appropriate
    underscore version of this function. Requires `n \ge 1`.

.. function:: void _acb_dirichlet_refine_hardy_z_zero(arb_t res, const arf_t a, const arf_t b, slong prec)

    Sets *res* to the unique zero of the Hardy Z-function in the
    interval `(a, b)`.

.. function:: void acb_dirichlet_hardy_z_zero(arb_t res, const fmpz_t n, slong prec)

    Sets *res* to the *n*-th zero of the Hardy Z-function, requiring `n \ge 1`.

.. function:: void acb_dirichlet_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, slong prec)

    Sets the entries of *res* to *len* consecutive zeros of the
    Hardy Z-function, beginning with the *n*-th zero. Requires positive *n*.

.. function:: void acb_dirichlet_zeta_zero(acb_t res, const fmpz_t n, slong prec)

    Sets *res* to the *n*-th nontrivial zero of `\zeta(s)`, requiring `n \ge 1`.

.. function:: void acb_dirichlet_zeta_zeros(acb_ptr res, const fmpz_t n, slong len, slong prec)

    Sets the entries of *res* to *len* consecutive nontrivial zeros of `\zeta(s)`
    beginning with the *n*-th zero. Requires positive *n*.

.. function:: void _acb_dirichlet_exact_zeta_nzeros(fmpz_t res, const arf_t t)

.. function:: void acb_dirichlet_zeta_nzeros(arb_t res, const arb_t t, slong prec)

    Compute the number of zeros (counted according to their multiplicities)
    of `\zeta(s)` in the region `0 < \operatorname{Im}(s) \le t`.

.. function:: void acb_dirichlet_backlund_s(arb_t res, const arb_t t, slong prec)

    Compute `S(t) = \frac{1}{\pi}\operatorname{arg}\zeta(\frac{1}{2} + it)`
    where the argument is defined by continuous variation of `s` in `\zeta(s)`
    starting at `s = 2`, then vertically to `s = 2 + it`, then horizontally
    to `s = \frac{1}{2} + it`. In particular `\operatorname{arg}` in this
    context is not the principal value of the argument, and it cannot be
    computed directly by :func:`acb_arg`. In practice `S(t)` is computed as
    `S(t) = N(t) - \frac{1}{\pi}\theta(t) - 1` where `N(t)` is
    :func:`acb_dirichlet_zeta_nzeros` and `\theta(t)` is
    :func:`acb_dirichlet_hardy_theta`.

.. function:: void acb_dirichlet_backlund_s_bound(mag_t res, const arb_t t)

    Compute an upper bound for `|S(t)|` quickly. Theorem 1
    and the bounds in (1.2) in [Tru2014]_ are used.

.. function:: void acb_dirichlet_zeta_nzeros_gram(fmpz_t res, const fmpz_t n)

    Compute `N(g_n)`. That is, compute the number of zeros (counted according
    to their multiplicities) of `\zeta(s)` in the region
    `0 < \operatorname{Im}(s) \le g_n` where `g_n` is the *n*-th Gram point.
    Requires `n \ge -1`.

.. function:: slong acb_dirichlet_backlund_s_gram(const fmpz_t n)

    Compute `S(g_n)` where `g_n` is the *n*-th Gram point. Requires `n \ge -1`.

Riemann zeta function zeros (Platt's method)
-------------------------------------------------------------------------------

The following functions related to the Riemann zeta function use the ideas
and formulas described by David J. Platt in [Pla2017]_.

.. function:: void acb_dirichlet_platt_scaled_lambda(arb_t res, const arb_t t, slong prec)

    Compute `\Lambda(t) e^{\pi t/4}` where

    .. math ::

        \Lambda(t) = \pi^{-\frac{it}{2}}
                         \Gamma\left(\frac{\frac{1}{2}+it}{2}\right)
                         \zeta\left(\frac{1}{2} + it\right)

    is defined in the beginning of section 3 of [Pla2017]_. As explained in
    [Pla2011]_ this function has the same zeros as `\zeta(1/2 + it)` and is
    real-valued by the functional equation, and the exponential factor is
    designed to counteract the decay of the gamma factor as `t` increases.

.. function:: void acb_dirichlet_platt_scaled_lambda_vec(arb_ptr res, const fmpz_t T, slong A, slong B, slong prec)

.. function:: void acb_dirichlet_platt_multieval(arb_ptr res, const fmpz_t T, slong A, slong B, const arb_t h, slong J, slong K, slong sigma, slong prec)

.. function:: void acb_dirichlet_platt_multieval_threaded(arb_ptr res, const fmpz_t T, slong A, slong B, const arb_t h, slong J, slong K, slong sigma, slong prec)

    Compute :func:`acb_dirichlet_platt_scaled_lambda` at `N=AB` points on a
    grid, following the notation of [Pla2017]_. The first point on the grid
    is `T - B/2` and the distance between grid points is `1/A`. The product
    `N=AB` must be an even integer. The multieval versions evaluate the
    function at all points on the grid simultaneously using discrete Fourier
    transforms, and they require the four additional tuning parameters
    *h*, *J*, *K*, and *sigma*. The *threaded* multieval version splits the
    computation over the number of threads returned by
    *flint_get_num_threads()*, while the default multieval version chooses
    whether to use multithreading automatically.

.. function:: void acb_dirichlet_platt_ws_interpolation(arb_t res, arf_t deriv, const arb_t t0, arb_srcptr p, const fmpz_t T, slong A, slong B, slong Ns_max, const arb_t H, slong sigma, slong prec)

    Compute :func:`acb_dirichlet_platt_scaled_lambda` at *t0* by
    Gaussian-windowed Whittaker-Shannon interpolation of points evaluated by
    :func:`acb_dirichlet_platt_scaled_lambda_vec`. The derivative is
    also approximated if the output parameter *deriv* is not *NULL*.
    *Ns_max* defines the maximum number of supporting points to be used in
    the interpolation on either side of *t0*. *H* is the standard deviation
    of the Gaussian window centered on *t0* to be applied before the
    interpolation. *sigma* is an odd positive integer tuning parameter
    `\sigma \in 2\mathbb{Z}_{>0}+1` used in computing error bounds.

.. function:: slong _acb_dirichlet_platt_local_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, const fmpz_t T, slong A, slong B, const arb_t h, slong J, slong K, slong sigma_grid, slong Ns_max, const arb_t H, slong sigma_interp, slong prec)
.. function:: slong acb_dirichlet_platt_local_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, slong prec)
.. function:: slong acb_dirichlet_platt_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, slong prec)

    Sets at most the first *len* entries of *res* to consecutive
    zeros of the Hardy Z-function starting with the *n*-th zero.
    The number of obtained consecutive zeros is returned. The first two
    function variants each make a single call to Platt's grid evaluation
    of the scaled Lambda function, whereas the third variant performs as many
    evaluations as necessary to obtain *len* consecutive zeros.
    The final several parameters of the underscored local variant have the same
    meanings as in the functions :func:`acb_dirichlet_platt_multieval`
    and :func:`acb_dirichlet_platt_ws_interpolation`. The non-underscored
    variants currently expect `10^4 \leq n \leq 10^{23}`. The user has the
    option of multi-threading through *flint_set_num_threads(numthreads)*.

.. function:: slong acb_dirichlet_platt_zeta_zeros(acb_ptr res, const fmpz_t n, slong len, slong prec)

    Sets at most the first *len* entries of *res* to consecutive
    zeros of the Riemann zeta function starting with the *n*-th zero.
    The number of obtained consecutive zeros is returned. It currently
    expects `10^4 \leq n \leq 10^{23}`. The user has the option of
    multi-threading through *flint_set_num_threads(numthreads)*.