acb_dirichlet.h

Initialise t for the roots of unity of order n = order, optimised for num evaluations at precision prec, and free it.

Sets z to e^{2\pi i n / \text{order}}.

Sets res to \chi(n), which is 0 if \gcd(n, q) > 1 and a root of unity otherwise.

Sets v to (\chi(0), \chi(1), \ldots, \chi(\text{nv} - 1)), sharing the powers of one root of unity.

Set res to the value \chi_q(m, n) of the Conrey pairing, that is, the value at n of the character with Conrey number m; the second form takes the two characters as input.

Set res to the Gauss sum G_q(\chi) = \sum_{x \bmod q} \chi(x) e^{2\pi i x/q}. The naive version sums the definition; factor writes the sum as a product of local Gauss sums by the Chinese remainder theorem; order2 assumes that \chi is real and primitive and returns i^\delta \sqrt{q}; theta assumes that \chi is primitive and obtains the sum from the functional equation of the theta series at t = 1. The default version reduces an imprimitive character to the primitive one inducing it and then chooses among the methods.

Set res to the Jacobi sum J_q(\chi_1, \chi_2) = \sum_{x \bmod q} \chi_1(x) \chi_2(1-x). The naive version sums the definition; factor writes it as a product of local Jacobi sums; gauss assumes that \chi_1 \chi_2 is primitive and uses J_q(\chi_1, \chi_2) G_q(\chi_1\chi_2) = G_q(\chi_1) G_q(\chi_2). The default version handles trivial and inverse pairs by closed formulas (and returns 0 for even q) and otherwise chooses among the methods. The ui version takes the Conrey numbers a and b.

Set res to \Theta_q(\chi, t) for real t > 0; the ui version takes the Conrey number a. The series converges quickly only for t not too small; for t < 1 the caller should use the functional equation above, which is not applied automatically (to avoid recomputing the Gauss sum).

Return a number of terms N such that the tail of the theta series of modulus q at t is less than 2^{-\text{prec}}. The first version uses a lower bound for t.

Set res to \sum_{n<\text{len}} n^p z^{a_n} x^{n^2}, where p is parity (0 or 1), z is the root of unity of z, and the exponents a_n are given (with a special value for \chi(n) = 0, as produced by dirichlet_chi_vec_order). The naive version sums the series as written; the smallorder version sums it in the quotient ring by the cyclotomic polynomial and evaluates at the root of unity only at the end, ignoring z except for its order.

Sets res to \sum_{k=1}^{\text{len}-1} a_k x^{k^2} (the entry a_0 is ignored), computing x^{k^2} by two multiplications per term.

Sets v to (x^{k^2})_{0 \le k < \text{nv}}.

Sets res[k], for 0 \le k < n, to t^{k/2} \sum_{m=1}^{\text{len}-1} \chi(m) m^\delta e^{-\pi m^2 t^{2k}/q} (half the theta series at the geometric sequence t^k, scaled), with the tail bound for the truncation after len terms. These are samples for computing the Mellin transform of the theta series.

Sets res to k^{-(s+x)} as a power series in x of length len. The flags integer and critical_line select optimisations for integer s and for \operatorname{Re}(s) = 1/2.

Sets z to \sum_{k=1}^n k^{-(s+x)} as a power series in x of length len. Powers of composite k = ij are computed as i^{-s} j^{-s} from a table, and the even terms are grouped as a polynomial in 2^{-(s+x)}. This needs about n / \log n powers, n/2 multiplications and storage for n/6 power series, so it is only faster than the naive sum for small len.

Sets z to \sum_{k=1}^n k^{-(s+x)} as a power series in x of length len, by partial sieving that adds the multiples of 5-smooth numbers into separate buckets. It computes about 4/15 of the powers, more than acb_dirichlet_powsum_sieved, but needs only logarithmic extra space and is faster for large len, since power series multiplications are replaced by additions.

Sets res to \zeta(s), choosing the algorithm from s and prec: exact values at integers (Bernoulli numbers for negative integers, arb_zeta_ui for positive ones), the Riemann–Siegel formula for |\operatorname{Im}(s)| \ge 24 \, p^{3/2} and |\operatorname{Re}(s)| \le 10 + p/10 at precision p, and Euler–Maclaurin summation otherwise, after the reflection formula for \operatorname{Re}(s) < 0. The result is indeterminate if s contains the pole s = 1.

Sets res (a vector of length len) to the Taylor coefficients of \zeta(s + x), or of the deflated function \zeta(s + x) - 1/(s + x - 1) if deflate is nonzero. Euler–Maclaurin summation is used (with the reflection formula for \operatorname{Re}(s) < 0), or the Riemann–Siegel formula for len \le 2 high on the critical line.

Sets res to an upper bound for |\zeta(s)|, computed quickly: in and slightly around the critical strip by formula (43.3) of [Rad1973], to the right by the value at \operatorname{Re}(s), and to the left by the functional equation.

Sets der1 and der2 to upper bounds for |\zeta'(s)| and |\zeta''(s)|, meant for the critical strip. The bounds are not tight.

Sets res to the Dirichlet eta function \eta(s) = \sum_{k=1}^{\infty} (-1)^{k+1} k^{-s} = (1-2^{1-s}) \zeta(s) (the alternating character is not a Dirichlet character). The removable singularity at s = 1, where \eta(1) = \log 2, is handled.

Sets res to the Riemann xi function \xi(s) = \frac{1}{2} s (s-1) \pi^{-s/2} \Gamma(\tfrac{s}{2}) \zeta(s), which satisfies \xi(1-s) = \xi(s).

Sets c to F^{(j)}(p) for 0 \le j < N, by power series division. This breaks down at p = \pm 1/2, which does not occur for exact s.

Updates d in place from d_j^{(k-1)} to d_j^{(k)} for 0 \le j \le \lfloor 3k/2 \rfloor + 1 (for k = 0, initialises it).

Sets err to a bound for the error term RS_K, by Theorem 4.2 of [Ari2011].

Sets res to \mathcal{R}(s) for s in the upper half-plane, with exactly K asymptotic terms if K > 0 and an automatic choice otherwise.

Sets res to \zeta(s) by the Riemann–Siegel formula, with exactly K asymptotic terms if K > 0 and an automatic choice otherwise.

Sets res to the first len Taylor coefficients of \zeta at s by the Riemann–Siegel formula. Only len \le 2 is supported; the derivative is computed with a finite difference.

Sets res to the Hurwitz zeta function \zeta(s, a) = \sum_{k=0}^{\infty} (k + a)^{-s}, analytically continued. The case a = 1 calls acb_dirichlet_zeta, and s = 0 and small negative integers s use Bernoulli polynomials; otherwise Euler–Maclaurin summation is used (_acb_poly_zeta_cpx_series). The result is indeterminate if s contains 1.

Precomputes a grid of N Taylor polynomials of length K for \zeta(s, a + A), for fast evaluation of \zeta(s, a) at a \in (0,1] with fixed s. This costs NK Hurwitz zeta values; each later evaluation costs 2K arithmetic operations plus A powers, with a truncation error O(1/(2AN)^K). With deflate nonzero, the deflated function (without the pole at s = 1) is used.

Choose A, K and N to minimise the cost of num_eval evaluations at s and precision prec, and (for the first function) initialise pre with them. If evaluating each value from scratch is estimated to be cheaper, A = K = N = 0.

Frees the precomputed data.

Sets res to a bound for the truncation error (not including rounding errors) of the evaluation of \zeta(s, a) with parameters A, K, N, for 0 < a \le 1.

Sets res to \zeta(s, p/q) using the precomputed data, for 0 < p/q \le 1. Throws an exception if p > q.

Set res to the Lerch transcendent \Phi(z,s,a) = \sum_{k=0}^{\infty} z^k (k+a)^{-s}, analytically continued to |z| \ge 1. The direct version sums the truncated series; the integral version uses a Hankel contour integral; the default version checks special cases that reduce to the Hurwitz zeta function or polylogarithms and otherwise chooses between the two.

Sets res to the generalised Stieltjes constant \gamma_n(a) for an integer n \ge 0 (an exception is thrown for n < 0), the coefficient in the Laurent series \zeta(s,a) = \frac{1}{s-1} + \sum_{n=0}^\infty \frac{(-1)^n}{n!} \gamma_n(a) (s-1)^n. With a = 1 these are the ordinary Stieltjes constants.

Set res to the root number \epsilon(\chi) of a primitive character. The theta version uses \epsilon(\chi) = \Theta_q(\chi, 1) / \overline{\Theta_q(\chi, 1)}; the default version uses the Gauss sum, \epsilon(\chi) = G_q(\chi)/(i^\delta \sqrt q), when q has several prime factors, and the theta series otherwise. The default version throws an exception if \chi is not primitive.

Sets res to L(s,\chi), choosing the algorithm: the Riemann zeta function for q = 1, the functional equation for primitive \chi and \operatorname{Re}(s) < -1/2 (or near s = 0 for even \chi), the Euler product for large \operatorname{Re}(s), and otherwise the decomposition into Hurwitz zeta functions with a precomputation shared by the \phi(q) values.

Sets res to L(s,\chi) = q^{-s}\sum_{k=1}^q \chi(k) \zeta(s, k/q). For s = 1 and non-principal \chi, the deflated Hurwitz zeta function is used to avoid the poles. If precomp is NULL, each Hurwitz zeta value is computed directly; otherwise the precomputation is used.

Set res to L(s,\chi) by the Euler product, which is efficient for large \operatorname{Re}(s); the result is finite only if \operatorname{Re}(s) \ge 2. The underscore version takes an integer s and a real character given by the list chi of its values at 0, \ldots, \text{mod}-1, and computes 1/L(s, \chi) if reciprocal is nonzero.

Set res to L(s,\chi) for rational s. The afe version uses the approximate functional equation, a rapidly convergent series of incomplete gamma functions (evaluated by binary splitting and the bit-burst algorithm), which is fastest at very high precision or large modulus; the default version uses it when this is estimated to be faster and otherwise calls acb_dirichlet_l.

Sets res (of length \phi(q)) to L(s,\chi) for all characters \chi modulo q, from the \phi(q) values \zeta(s, k/q) and one discrete Fourier transform. The values are ordered by lexicographically ordered Conrey logarithms, as in acb_dirichlet_dft_index. If precomp is NULL, each Hurwitz zeta value is computed directly.

Sets res to the Taylor coefficients L(s), L'(s), \ldots, L^{(\text{len}-1)}(s)/(\text{len}-1)!. With deflate nonzero and principal \chi, the expansion of L(s,\chi) - \frac{\phi(q)}{q (s-1)} is computed instead, which at s = 1 is the regular part of the Laurent expansion; for non-principal \chi, deflate has no effect.

Set res to L(s(x),\chi) truncated to length len, for a power series s(x), with deflate as in acb_dirichlet_l_jet.

Sets w to the DFT of v, both vectors of length \phi(q) indexed by the lexicographically ordered Conrey logarithms (the order of dirichlet_char_next and dirichlet_char_index).

Sets w to the DFT of v, both vectors of length q indexed by Conrey numbers. The entries at indices not coprime to q are ignored.

Sets res to the first len Taylor coefficients of \theta(t), using the log-gamma function.

Sets res to the first len Taylor coefficients of Z(t), computed as e^{i\theta(t)} L(1/2 + it, \chi) with acb_dirichlet_l_jet. For the Riemann zeta function, this uses the Riemann–Siegel formula when t is large and len \le 2. If t is real, the imaginary parts are set to exactly zero.

Set res to \theta(t(x)), respectively Z(t(x)), truncated to length len, for a power series t(x) (given as s).

Sets res to the n-th Gram point g_n, the unique solution in [7, \infty) of \theta(g_n) = \pi n, for n \ge -1. Only the Riemann zeta function is supported: G and chi must be NULL; otherwise, or for n < -1, the result is indeterminate.

Returns a number B of consecutive good Gram blocks that suffices for Turing's method [Tur1953] with the bounds of [Leh1970], [Bre1979] and [Tru2011]: if 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.

Set res to N(t). For exact t the result is an exact integer; for an inexact ball, it is the ball containing N(t_1), \ldots, N(t_2) where [t_1, t_2] encloses t. The underscore version takes an exact t and sets an integer.

Sets res to N(g_n), for n \ge -1.

Sets res to S(t) = \frac{1}{\pi}\operatorname{arg}\zeta(\frac{1}{2} + it), where the argument is defined by continuous variation from s = 2 to 2 + it and then to \frac12 + it, not the principal value. It is computed as N(t) - \theta(t)/\pi - 1.

Sets res to an upper bound for |S(t)|, by Theorem 1 and (1.2) of [Tru2014].

Returns S(g_n), for n \ge -1.

Set res to \gamma_n, respectively to the len consecutive values \gamma_n, \ldots, \gamma_{n+\text{len}-1}, accurate to about prec bits. Throws an exception if n < 1. The zeros are isolated rigorously and computed by the Riemann–Siegel method, or by Platt's method at large height when that is expected to be faster.

Set res to \rho_n = 1/2 + i\gamma_n, respectively to len consecutive zeros starting with \rho_n. Requires n \ge 1.

The same as acb_dirichlet_hardy_z_zeros, always with the Riemann–Siegel method: isolation with Gram points and Turing's method, then refinement with _acb_dirichlet_refine_hardy_z_zero.

Returns 0 if acb_dirichlet_hardy_z_zeros uses the Riemann–Siegel method for len zeros from the n-th at target precision prec (or at the default precision 64 + \lceil \log_2 n \rceil if prec is not positive), and otherwise the working precision of the large height method. The choice minimises a cost model measured with 2 threads.

Set (a, b) to an interval that contains \gamma_n and no other zero of Z. The gram version uses Gram's law and requires 1 \le n \le 126; the rosser version uses Rosser's rule and requires 1 \le n \le 13999526; the turing version uses Turing's method and requires n \ge 2; the default version chooses among them and requires n \ge 1.

Sets res to Z(t), increasing the precision *pprec as needed until the sign of Z(t) is determined, and returns the sign.

Sets res to the unique zero of Z in (a, b), to prec bits.

Sets res to the zero of Z in the ball z, which must contain exactly one zero (for example a zero computed by Platt's method), to prec bits: z rounded if it is accurate enough, otherwise the last stage of the refinement above if z is accurate to about half the precision, and the full refinement otherwise.

Sets res to \Lambda(t) e^{\pi t/4}, where \Lambda(t) = \pi^{-it/2} \Gamma(\frac{1/2+it}{2}) \zeta(\frac12 + it) is defined at the beginning of §3 of [Pla2017]. By the functional equation it is real for real t, it has the same zeros as \zeta(1/2 + it) [Pla2011], and the exponential factor cancels the decay of the gamma factor.

Sets res to the values of acb_dirichlet_platt_scaled_lambda at the N = AB grid points T - B/2 + k/A, 0 \le k < N, by evaluating each point separately. N must be even.

Set out to the values of acb_dirichlet_platt_scaled_lambda on the same grid as acb_dirichlet_platt_scaled_lambda_vec, all at once with discrete Fourier transforms, in the notation of [Pla2017]: h is the width of the Gaussian, J the number of terms of the Dirichlet series, K the number of Taylor terms and sigma an odd integer used in the error bounds. The threaded version uses flint_get_num_threads() threads; the default version decides whether to use threads.

Sets res to acb_dirichlet_platt_scaled_lambda at t_0, by Gaussian-windowed Whittaker–Shannon interpolation of the grid values p. If deriv is not NULL, it is set to an approximation of the derivative. Ns_max is the maximum number of grid points used on either side of t_0, H is the standard deviation of the Gaussian window, and sigma is an odd positive integer used in the error bounds.

Set at most len entries of res to consecutive zeros \gamma_n, \ldots of the Hardy Z-function and return the number found. The local versions make one grid evaluation; acb_dirichlet_platt_hardy_z_zeros repeats them until len zeros are found. The parameters of the underscore version have the meanings of acb_dirichlet_platt_multieval and acb_dirichlet_platt_ws_interpolation.

Like the local functions above, but set res to isolating intervals of the zeros instead of refining them.

Sets at most len entries of res to consecutive zeros \rho_n, \ldots of \zeta(s) by Platt's method, and returns the number found.

Sets res to the secondary zeta function Z(s) = \sum_{n=1}^{\infty} \gamma_n^{-s}, where \rho_n = 1/2 + i\gamma_n runs over the zeros of \zeta(s) with positive imaginary part. The Riemann hypothesis is not assumed. The algorithm follows [Ari2020], with the decomposition Z(s) = A(s) - P(s) + E(s) - S(s) depending on a parameter a > 0, an improved tail bound for P, and rigorous bounds for the tails of A and S that replace the heuristic bounds of [Ari2020].

Sets xt to \pi t^2 / q, the exponent of the theta series.

Sets res to \sum_{n=1}^{\text{len}-1} \chi(n) n^\delta x^{n^2} for x = xt, using acb_dirichlet_qseries_arb.

Sets res to an upper bound for \sum_{k > n} k e^{-a k^2}.

Sets res to the polynomial with the len coefficients v evaluated at z, by the baby-step giant-step method.

Returns the local factor at p (with p^e \| q and pe = p^e) of the Jacobi sum J_q(1, \chi) for a character of conductor cond. The value is a signed integer cast to ulong.

Sets res to \beta(t) = 1/6 + \log\log t / \log t, the exponent in the bound |\zeta(1/2 + it)| = O(t^{\beta(t)}) used in [Pla2017].

Set out to the error bounds of Lemma 3.2 and Lemmas A.5, A.7, A.9, A.11, B.1 and B.2 of [Pla2017], which bound the errors of the grid evaluation in acb_dirichlet_platt_multieval (truncation of the Dirichlet series, of the Taylor expansion, aliasing of the discrete Fourier transform, and so on). Lemma A.5 requires B > h\sqrt{k}; when a lemma's hypotheses fail, the bound is infinite.

Set res to the bound C(\sigma, t_0, h, k) of [Pla2017] for the k-th derivative of the scaled function in the strip used for the error analysis, with or without precomputed constants that depend only on \sigma, h and k.

Set res to the bound I of [Pla2017] for the error of the Whittaker–Shannon interpolation caused by the Gaussian window and the sampling rate A, with or without precomputed constants.

The interpolation acb_dirichlet_platt_ws_interpolation with the constants of both bounds precomputed, for many evaluations with the same A, H and \sigma.

Sets res to the bound C_3 of [Pla2017] for the error of truncating the interpolation sum to N_s points on either side.