Series with Binomial-like Coefﬁcients for the Investigation of Fractal Structures Associated with the Riemann Zeta Function

: The paper continues the study of efficient algorithms for the computation of zeta functions over the complex plane. We aim to apply the modifications of algorithms to the investigation of underlying fractal structures associated with the Riemann zeta function. We discuss the computational complexity and numerical aspects of the implemented algorithms based on series with binomial-like coefficients.


Introduction
In this paper, we continue the study of efficient algorithms for computation of the Riemann zeta function over the complex plane, introduced by Borwein [1] and extended by Belovas et al., (see [2,3] and references therein). Šleževičienė [4], Vepštas [5], and Coffey [6] applied this methodology for the computation of Dirichlet L-functions, Hurwitz zeta function, and polylogarithm. Belovas et al. obtained limit theorems, which allowed the introduction of asymptotic approximations for the coefficients of the series of the algorithms. A preliminary presentation of computational aspects of the approach has been presented in [3]. Theoretical aspects of the approach (as well as more subtle proofs of the limit theorems) have been discussed in [7].
Fractal geography of the Riemann zeta function (and other zeta functions) was addressed by King [8]. Woon [9] and Tingen [10] computed Julia and Mandelbrot sets of the Riemann zeta function and Hurwitz zeta function, respectively, and studied the properties of these fractals. Recently Blankers et al. [11] investigated the analogs of Julia and Mandelbrot sets for dynamical systems over the hyperbolic numbers. In the present study, we enhance algorithms for the calculation of the Riemann zeta function, proposed in [2,3]. We specify the convergence rate to the limiting distribution for the coefficients of the series, identify the error term, and discuss computational complexity. The algorithms are compared against the recently proposed Zetafast algorithm [12] and are applied for the investigation of underlying fractal structures associated with the Riemann zeta function.
The paper is organized as follows. The first part is the introduction. In Section 2, we describe algorithms and present theoretical results. Section 3 is devoted to the visual investigation of the underlying fractal background of the Riemann zeta function. Pseudocodes of the algorithms for the computation and the visualization are given in Section 4. Sections 5 and 6 are devoted to presenting the results and conclusions, respectively.
Throughout this paper, U × V stands for the Cartesian product of sets U and V. We denote by Φ(x) the cumulative distribution function of the standard normal distribution, and by Γ(s) we denote the gamma function. Next, x and x stand for the floor function and the ceiling functions, respectively. All limits in the paper, unless specified, are taken as n → ∞.

MB-and BLC-Algorithms for the Computation of the Riemann Zeta Function
Let s = σ + it be a complex variable. The Riemann zeta-function is defined on the half-plane σ > 1 by the ordinary Dirichlet series or the Euler product formula, and by analytic continuation for other complex values. The Riemann zeta function is a meromorphic function (holomorphic on the whole complex plane except for a simple pole at s = 1 with residue 1). The Riemann zeta function satisfies the functional equation implying that ζ(s) has simple zeros at s = −2n, n ∈ N, known as the trivial zeros.
Other zeroes are called nontrivial. The famous Riemann hypothesis states that all the nontrivial zeros lie on the critical line s = 1/2 + it. The hypothesis is closely related to the distribution of prime numbers, implying the best possible error term in the prime number theorem, Here π(x) is the prime-counting function, i.e., the number of primes less than or equal to x , x ∈ R. A summary of the literature covering the problems related to the Riemann zeta function and its applications is presented in [13,14] and the references therein.

MB-Algorithm
In [3] Belovas et al., proposed a modification of Borwein's efficient algorithm (MB-algorithm) for the computation of the Riemann zeta function [1]. The algorithm applies to complex numbers with σ 1/2 and arbitrary t. The Riemann zeta function is represented by the alternating series Here (case j = 1 in ψ (j) n,k corresponds MB-series), by Theorem 1 from [3], we have while The algorithm is nearly optimal in the sense that there is no sequence of n-term exponential polynomials that converge to the Riemann zeta function much faster than of the algorithm (see Theorem 3.1 in [1]).
Here I x (a, b) stands for the regularized incomplete beta function, The error terms γ (j) n (s) of these methods are discussed in the following subsection.

Error Terms and Computational Complexity
First we formulate an auxiliary lemma, aiming to investigate the behaviour of the series in the neighbourhoods of critical points τ k , Note that in (1) the denominator 1 − 2 1−s = 0 if and only if s = 2πk/log 2, k ∈ Z and s = 1.
The error term and the computational complexity are closely linked to the problem of the selection of the minimal number of terms in the series (1). Let us formulate the following theorem. Theorem 1. Let σ 1/2, t 0, ε > 0 and |s − τ k | ε, then (i) the error term of the series (1) is (ii) the series (1) to compute the Riemann zeta-function with d decimal digits of accuracy, require a number of terms with coefficients of expressions (9) and (10) presented in Table 1. Table 1. Coefficients of expressions (9) and (10).
Considering the function I(σ), we have By a product representation of the gamma function (cf. 8.326.1 in [15]), The product is decreasing by σ, hence (cf. 8.332.2 in [15]), Hence, In view of (14), to compute the Riemann zeta-function with d decimal digits of accuracy, the approach requires a number n of terms not less than 2 0 . Let |s − τ k | ε and |σ − 1| ε. By applying the maximum modulus principle and Lemma 1, we receive thus concluding the proof. The deduction for BLC-series is analogical.

Corollary 1.
Under the conditions of Theorem 1, for ε = 10 −m , m ∈ N, the series (1) to compute the Riemann zeta-function with d decimal digits of accuracy, requires the number of terms Proof of Corollary 1. The result (18) follows immediately, if we notice that for ε → 0 we have

NA-Modifications of MB-and BLC-Algorithms
Limit theorems for coefficients of MBand BLC-series enable us to derive a normal approximation for coefficients ψ (j) n,k (cf. (24) in [3]). We can formulate the following theorem.
Coefficients µ (j) n and σ (j) n are presented in Table 2.
Proof of Theorem 2. Let us start with MB-series coefficients. Suppose A n is an integral random variable with the probability mass function Here (cf. (1) in [3]) Thus, Let F n (x) be the cumulative distribution function of the random variable A n (20), then (cf. Theorem 3 in [7]) Note that the cumulative distribution function Denoting k = σ n x + µ n and taking into account (22) and (23), we obtain The first part of the theorem follows. Similar result for BLC-coefficients ψ (2) n,k has been proven in [2].
Theorem 2 allows us to choose the number of terms n (j) for the series (1), for n large enough. Here for fixed σ and d. The refined version of N A-modification based methodology is summarized in Section 4.

Empirical Insights for NA-Modifications
While performing practical computations using N A-algorithms, we have noticed that the values produced were significantly more accurate than otherwise implied by d in the analytic estimate (10). In order to increase the performance and to have a clear course for future theoretical refinements, we propose empirical formulae for the minimum number of terms in the series (1) to compute the Riemann zeta-function with d decimal digits of accuracy.
Kuzma has proposed the following empirically-based estimate for the number of terms for the BLC-series (d = 6), In the present section, we offer an improvement to this estimate. Figure 1 displays the minimum n required to calculate the Riemann zeta function with d = 6 digits of accuracy using N Aand BLC-algorithms at σ = 1/2, t ∈ [1000, 1050] (the blue curve). The curves have clearly visible periodic peaks (marked by red vertical lines). The peaks have a period of λ = 2π/ log 2, which correspond τ k special points of Theorem 1. Since we are interested in the upper bound of this empirical curve, for the following calculations we use the points t = λk, k ∈ N.
derived for d ∈ [1, 10] using the points (σ, t) ∈ 1/2 × (0, 10,000). Each graph represents a fitted curve for a different d value.  Here we can clearly see that a (1) has no correlation with d while b (1) and with c (1) = xd + y we obtain the following coefficients for (26) (see Table 3):

Methods of the Visualization
In this study we employ two methods to reveal the Riemann zeta function underlying nature. The first heuristic method (FH-method) calculates RGB colors of the graph of the Riemann zeta function, using a composition of special functions. Suppose we have a function f : (R + , C) → N 0 : Now we can define functions f 1 , f 2 , f 3 : Next, we calculate (R, G, B) colors of each pixel of the graph of the Riemann zeta function using polynomial functions of f k (see Table 4): The second approach (second fractal heuristic (SFH) method) is based on the application of the Mandelbrot set to the visualization of the Riemann zeta function. Suppose we aim to visualize ζ(σ + it) for (σ, t) ∈ (σ 1 , σ 2 ) × (t 1 , t 2 ). First, we introduce the logtransformation for each point (x, y) of the graph, thus obtaining the set Q = (x min , y min ) × (x max , y max ). Here Next, we linearly transform Q into the subset S of the complex plane, We take S = (−2, 0.47) × (−1.12i, 1.12i), where the Mandelbrot set is defined. Then we use an algorithm to generate the Mandelbrot set, setting the start position at z 0 = 0 and z * = (x * , y * ): Suppose that k ∈ N, k v max indicates the number of iterations (31), required to ascertain that z * does not belong to the Mandelbrot set, with |z k+1 | 2 and k < v max .
For k = v max , it is unclear if z * does not belong to the Mandelbrot set. Now let k 0 = 50k . We calculate RGB color for the z * point by the following rule:

Visual Investigations
The first visualization (see Figure 4) reveals the underlying structures in the "center" S 1 ⊂ C of the Riemann zeta function, received by two different methods (the color visualization and the fractal visualization). Here S 1 = (−20, 8) × (−14, 14). Figure 4a is obtained using FH-method with color parameters η 1 = 100 and η 2 = η 3 = 8. The color transform g (1) k is linear (see Table 4). Figure 4b is obtained using SFH-method. Note small bright fractal feature on the right-hand side, calling for in-depth investigation. Figure 5 presents zoom-in frames of S 2 region for the Riemann zeta function. Here S 2 = (−5, 6) × (β, α + β), with four shifted in β intervals (see Table 5 for the ranges). The frames were received using FH-method with color parameters η 1 = 100 and η 2 = η 3 = 8. The color transform g (1) k is linear (see Table 4). Note nontrivial zeros of the Riemann zeta function (blue disks, marked with arrows in Figure 5a   , received by SH and SFH methods. Note small fractal feature on the right-hand side of (b).
(a) (b) (c) (d) Figure 5. FH-based zoomed-in frames of the Riemann zeta function (see Table 5 for the ranges). Note nontrivial zeros of the Riemann zeta function (blue disks, marked with arrows in Figure 5a,b). Table 5. Ranges of the sets of Figure 5: (σ, t) ∈ S 2 , α = 50.   Figure 6. Fractal features of the Riemann zeta function in the pole area (see Table 6 for the ranges). (a) gives zoomed-in image of the feature observed in Figure 4b. (b) represents the next magnification step (red rectangle). Table 6. Ranges of the sets of Figure 6, (σ, t) ∈ S 3 .    Figure 6b. The next five frames (each of them corresponds to a colored rectangle in Figure 7a) uncover some aesthetically pleasing features of fractal structures associated with the Riemann zeta function. Note snowflake-shaped fractals in Figure 7c, as well as pinwheel-shaped ones in Figure 7d,e, resembling discs of spiral galaxies. Clockwise spinning Figure 7e reminds us of the grand design spiral galaxy NGC 4254 in Coma Berenices. Counter-clockwise rotating Figure 7d resembles the Pinwheel Galaxy NGC 5457 in Ursa Major. Invariant features of fractal geometry generated from images provide a good set of descriptive values for the recognition of regions and objects, e.g., fractal signatures of galaxies are examined with the aim of classifying them (cf. [16]). Figure 7 is received by SFH-method. The ranges of the sets are given in Table 7.   Table 7. Table 7. Ranges of the frames of Figure 6.
Color parameters are given in Table 8. Detailed color parameters are given in Table 8. Table 8. Color parameters of Figure 8.

Computation and Visualization Algorithms
This section gives pseudocodes of the algorithms described in Sections 2 and 3.

Computation Algorithms
The first algorithm outlines MBand BLC-approaches (cf. Theorem 1 and Corollary 1) with the corresponding empirical modifications (26) for the calculation of multiple values of the Riemann zeta functions while t is fixed.
The second algorithm outlines N A-modifications of MBand BLC-methods. These approaches are more suitable for the calculation of specific values of the Riemann zeta function.
Results of numerical experiments with Algorithms 1 and 2 are presented in Section 5.
Algorithm 1 This algorithm will return multiple values of the Riemann zeta function for fixed t and array {σ r }. Note that L k = log k stand for precalculated logarithms.

Visualization Algorithms
The third algorithm (Algorithm 3), corresponding the first heuristic method (FH-method), calculates RGB colors of the graph of the Riemann zeta function, using a composition of special functions.

Algorithm 4
This algorithm will return fractalized image of Riemann zeta function for (σ, t) ∈ (σ min , σ max ) × (t min , t max ). Here m stands for max iterations to get more precise fractal image, w-width in pixels of output image img. The output image utilizes yellow-black-blue color palette. 1: procedure SFH(σ min , σ max , t min , t max : real numbers; w, m : natural numbers) 2:

Numerical Experiments
We have performed numerical experiments with seven methods and modifications listed in Table 9.

First Numerical Experiment
The first numerical experiment deals with normal approximation-based modifications (cf. Algorithm 2). Using N AMB (j = 3), N ABLC (j = 4) and Zetafast (j = 7) methods we generate sequences of values of the Riemann zeta function {ζ where s k p stand for critical points (5) with k p = 2 p+6 , 1 p 3, and ρ 1 = 10 −1 . Thus we obtain 9 sequences overall (3 algorithms × 3 sets of arguments). Using Zetafast algorithm as a benchmark we calculate the accuracy δ (j) p and the relative performance θ where τ (j) p is the processing time of jth sequence {ζ (j) l,p }, 1 l N, for fixed p. The results of the first numerical experiment are presented in Table 10.

Second Numerical Experiment
The second numerical experiment aims to verify the accuracy of the algorithms on fixed horizontal lines, close to critical points. Using MB (j = 1) and BLC (j = 2) methods, their empirical modifications (j = 5 and j = 6) and Zetafast method (j = 7), we generate (cf. Algorithm 1) sequences of values of the Riemann zeta function {ζ Thus we obtain 15 sequences overall (5 algorithms × 3 sets of arguments). Using the Zetafast algorithm as a benchmark we calculate the accuracy δ (j) p and the relative performance θ (j) p (cf. (33)). The results of the second numerical experiment are presented in Table 11.
The numerical experiments have been performed on Intel ® Core TM i7-8750H 2.2 GHz (boosted to 4.0 GHz) processor with 16 GB DDR4 RAM. The code has been compiled with g++ 11.2.0 compiler using O3 optimization. C++ Boost library has been used for the implementation of the incomplete beta function for BLC-algorithm. Table 11. Results of the second numerical experiment: accuracy δ (j) p and relative performance θ (j) p on fixed lines t p , for d = 6, m = 1. The last line shows the performance of ZF-algorithm (sec).
We have refined the error terms and the expressions for the minimal number of terms in MBand BLC-series of efficient algorithms for the computation of the Riemann zeta function, taking into account the behavior of the series in the neighborhoods of critical points. Theorem 1 shows that MB-based algorithms converge faster than BLC-based algorithms.
Indeed, the MB-coefficient of the error term Θ  (9)). However, BLC-approach has its advantages that might be useful in analytical research (cf. (4)). Note that this deficiency of the MB-algorithm is solved by the introduction of N A-modification (19).
The results of the numerical experiments (see Tables 10 and 11) show that MB and BLC methods, along with their normal and empirical modifications, allow fast and accurate calculations of the Riemann zeta function for large values of argument t. The results demonstrate that the introduced modifications accelerate computations of the Riemann zeta function, compared to Zetafast method. These versions of algorithms are well-suited for distributed computations and grid computing.

Findings of Visual Investigations of Fractal Structures, Associated with the Riemann Zeta Function
The illustrations obtained using FH-method clearly show the arrangement of trivial and non-trivial zeros of the Riemann zeta function in the complex plane (see Figure 5a,b). In addition to these points, we can also see dark 2D curves that satisfy the conditions (ζ(σ + it)) = 0 and (ζ(σ + it)) = 0 (see Figure 4a). The SFH-method distributes deformed copies of the Mandelbrot set in the complex plane, thus relating the values of the Riemann zeta function to the fractal structure. This allows for a visual assessment of essential changes in the Riemann zeta function values. Next, SFH-approach reveals notable symmetric fractals characterizing the neighborhood of the pole of the Riemann zeta function (see Figures 6 and 7).

Future Research Directions
Numerical experiments with empirical formulas indicate that the theoretical selection of the number of terms of the series n can be reduced. Next, the accuracy of the normal approximationbased modifications of MB and BLC algorithms might be refined by employing the theory of large deviations. The figures presented in this work reveal areas of the complex plane where the modulus of the Riemann zeta function exhibits very volatile values. This allows us to investigate the complex plane regions of ζ(s) = ζ(s), thus enabling us to locate non-trivial zeros' positions visually. In future works, these visual instruments could be refined.