Sieved power sums with dfloat balls

These functions compute the partial sums

that make up the main sum of the Riemann–Siegel formula for $\zeta(s)$ at large height $t = \operatorname{Im}(s)$. gr_powsum_sieved is a generic implementation over three rings given at run time, and dfloat_powsum_sieved instantiates it with dfloat balls and plain expansions. acb_dirichlet_zeta_rs_r calls dfloat_powsum_sieved for $N \ge 16$ and falls back to acb_dirichlet_powsum_sieved when it returns 0; this makes the Riemann–Siegel evaluations used to isolate zeros 17 to 18 times faster at the heights of the $10^{15}$-th to $10^{17}$-th zeros.

The algorithm#

The algorithm is that of acb_dirichlet_powsum_sieved. The function $f(n) = n^{-s}$ is completely multiplicative, so the term of a composite $n$ is the product $f(d) f(n/d)$ with $d$ the smallest prime factor of $n$, taken from a table of $f(k)$ for odd $k \le N/3$; only the prime terms need logarithms, sines, cosines and exponentials. The even $n$ are handled by Horner's rule in $2^{-s}$ over the dyadic ranges of the odd part.

The implementation rearranges this so that all the bulk work goes through vector operations of the generic rings:

• the odd numbers are processed in blocks of 4096 with a segmented sieve that carries, for each small prime, its next multiple and the cofactor of that multiple, so that no divisions are needed;
• the primes of a block are evaluated together: _gr_vec_log in the phase ring $P$, the phase $(t/2\pi)\log p$ reduced by subtracting a nearby integer exactly, rounded to the term ring $T$ with _gr_vec_set_other, then _gr_vec_sin_cos, _gr_vec_exp and _gr_vec_mul; on the critical line $\sigma = 1/2$ the magnitudes $p^{-1/2}$ come from _gr_vec_rsqrt instead of the exponential;
• the composite terms of a block are gathered from the table with _gr_vec_gather (real and imaginary parts interleaved, so that a lookup touches one cache line), multiplied with four _gr_vec_mul, and scattered back with _gr_vec_scatter;
• the terms of each dyadic range are summed pairwise with _gr_vec_add.

The phase ring $P$ needs about $\log_2(t \log N)$ more bits than the term ring $T$, because the phase $t \log p$ is reduced modulo $2\pi$.

Bounded memory: the two-phase version#

The table of $f(k)$ for $k \le N/3$ takes $8N/3$ bytes per double of the table ring, which limits the plain algorithm. With a table bound $N^{2/3} \le M < N/3$, the odd $k \le M$ are processed as above, and among the odd $c > M$ only the roots are computed: the primes and the composites with cofactor $c/d \le M$. Every $n \le N$ whose odd part exceeds $M$ is uniquely $n = mc$ with $c$ a root and $m \le N/c < N/M$ an integer whose odd part has no prime factor above $d$ (the chain of divisions of $n$ by its smallest prime factor reaches $c$ as its last element above $M$). Hence

where $G(z, p)$ is the sum of $f(m)$ over the $m \le z$ whose odd part has no prime factor above $p$, and $d_c = \infty$ for a prime $c$. With $Z = \lfloor N/(M+1) \rfloor \le N^{1/3}$, the values of $G$ form a table of about $Z^2/(2\log Z)$ entries. The second phase gathers $G(\lfloor N/c \rfloor, d_c)$ for each root and needs one complex product per root; all other composites (all but a fraction of about $1.12/\log(N/M)$) are never computed.

The default table bound (M = 0) is the full table for $N < 2048$ and $\lceil N^{2/3} \rceil$ otherwise. Besides bounding the memory by $O(N^{2/3})$, this is faster, since the table fits in the cache and fewer products are needed: at $N = 10^8$, with $(T, P) = (1, 3)$, the time per term drops from 24.6 to 15.8 ns and the memory from 262 MB to 11 MB.

The second phase, all but $O(N^{2/3})$ of the work, is split into chunks of at least $2^{18}$ odd numbers, each with its own sieve state, and run with flint_parallel_do on flint_get_num_threads threads. The partial sums are added in the order of the chunks, so the result does not depend on the number of threads.

Error bounds#

With F = NULL, every term is a ball and the result is rigorous by construction. With a plain ring $F$, only the prime terms are balls, and the error of the composite terms and block sums is bounded a priori and added per dyadic range. With $\varepsilon$ the largest relative radius of a prime term (obtained with _gr_vec_get_interval_mid_rad) and $\mu = \sqrt{2}(\delta_m + \delta_a + \delta_m \delta_a)$ the relative error of a plain complex product, a term $f(n)$ with $\Omega(n)$ prime factors has a relative error of at most $(1 + \varepsilon)^{\Omega}(1 + \mu)^{\Omega - 1} - 1$, where $\Omega \le \log_3 N$; the pairwise sums of at most 4096 terms add $(1 + \delta_a)^{12} - 1$ times the sum of the magnitudes, which is bounded by an integral. In the two-phase version, the sum over the roots of $g(c) = c^{-\sigma}$ times the sum of $g(m)$ over the corresponding $m$ is exactly the sum of $g(n)$ over the $n \le N$ with odd part above $M$, so the bound is as tight as in the first phase.

The dfloat instantiation uses $\delta_m = \delta_a = u$ for d1 and $\delta_m = 8u^2$, $\delta_a = 3u^2$ for d2, $u = 2^{-53}$, from the error analysis of the kernels (Double-word and multi-double arithmetic); the phase context is strong (canonical), so that rounding the reduced phases to $T$ is a truncation.

Choosing the precisions#

dfloat_powsum_sieved chooses $T$ so that the terms have $b = \mathit{prec} + \log_2(N)/2 + 8$ bits (the accumulated error grows like $\sqrt{N} u^T$, so that the result has about prec bits after the binary point), and $P$ so that the phase $t \log N$ leaves $b + 6$ bits after the binary point. With $P = 4$, a phase short by up to 8 bits is accepted. It returns 0 if no pair $(T, P)$ in $\{(1, 2), (1, 3), (2, 3), (2, 4), (3, 4)\}$ fits, if $|\operatorname{Re}(s)| > 2^{10}$, $|t| \ge 10^{60}$ or $N \ge 2^{50}$, if the tables would take more than DFLOAT_POWSUM_MAX_BYTES, or if dfloat_is_supported returns 0. The plain composites are used for $T \le 2$; _dfloat_powsum_sieved_ball keeps every term a ball (and is the only variant for $T = 3$): its radii are about 10 times smaller and it is 20 to 30% slower.

Performance#

At the height of the $10^{15}$-th zero ($t \approx 2.1 \cdot 10^{14}$, $N = 5\,760\,732$, one core), acb_dirichlet_powsum_sieved at 179 bits takes 2.4 s, while the full-table dfloat version takes 0.13 s with $(1, 2)$ (radius $2^{-30}$), 0.27 s with $(2, 3)$ (radius $2^{-79}$) and 0.60 s with $(3, 4)$ (radius $2^{-133}$); the default two-phase version is another 25 to 30% faster. About half of the time goes to the prime terms, a quarter to the composites (dominated by cache misses of the table lookups), and the rest to the sieve and the sums. At low heights the fixed setup cost of a few microseconds makes the arb_t version faster for $N < 16$, which is why acb_dirichlet_zeta_rs_r uses the dfloat version only from $N = 16$ on.

@misc{flint-doc-dfloat-powsum,
  author       = {{The FLINT team}},
  title        = {{Sieved power sums with dfloat balls (dfloat.h)}},
  howpublished = {FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev},
  year         = {2026},
  url          = {https://flintlib.org/doc/dfloat/powsum.html},
  note         = {Generated 2026-11-02}
}