The modified Borwein’s algorithm (MB-algorithm) applies to complex numbers with ℜ s ⩾ 1 / 2 and arbitrary ℑ s. In this modification (see Belovas et al., 2022), the Riemann zeta function is represented by the alternating series (1) ζ ( s ) = 1 1 − 2 1 − s ∑ k = 0 n − 1 ( − 1 ) k ψ n , k ( k + 1 ) s + R n ( s ) . Here ψ n , k = 1 − H k H n and l max = arg max 0 ⩽ k ⩽ n ( n + k − 1 ) ! 4 k ( n − k ) ! ( 2 k ) ! , n ∈ N , 0 ⩽ k ⩽ n , while H l = H l − 1 + exp ( T l − T l max + ( l − l max ) log 4 ) , H 0 = exp ( T 0 − T l max − l max log 4 ) , T l = T l − 1 + log ( n − l + 1 ) ( n + l − 1 ) ( 2 l − 1 ) ( 2 l ) , T 0 = − log n , 1 ⩽ l ⩽ n .

The error term and the computational complexity are closely linked to the problem of selecting the minimal number of terms in the series (1). Theorem 1 in Belovas et al. (2022) states that the series (1) to compute the Riemann zeta-function with d decimal digits of accuracy require a number of terms (2) n = ⌈ π t 2 log ( 3 + 8 ) + log 10 log ( 3 + 8 ) d + log 2 1 − 2 − ε log ( 3 + 8 ) ⌉ , here | s − s k | ⩾ ε > 0 and s k = 1 + 2 i π k / log 2.

The limit theorem for the coefficients of MB-series enables us to derive the following normal approximation for the coefficients ψ n , k of the series (1) (cf. Belovas et al., 2022) ψ n , k = 1 − Φ ( k − μ n σ n ) + O ( 1 n ) , μ n = n 2 , σ n = n 32 4 . Here Φ ( x ) is the cumulative distribution function of the standard normal distribution. This result allows us to select the number of terms k 0 in the series (1), (3) k 0 = ⌈ μ n + z σ n ⌉ , z = Φ − 1 ( 1 − 10 − d ) , n > 10 3 . Note that for fixed ℜ s and d, we have k 0 ∼ 0.63 ℑ s.

Kuzma (Belovas et al., 2022) proposed the following empirical (EMB) formula for the minimal number of terms n in the series (1) for calculations with d decimal digits of accuracy (4) n = ⌈ a t + b t + c ⌉ , a = 0.451 , b = 1.407 d − 0.245 , c = 0.371 d + 0.195. The empirical formula (4) has been verified for 0 ⩽ t ⩽ 10 4 .

Associating polynomials p n ( x ) in Proposition 1 with Chebyshev polynomials T n ( x ) = ( − 1 ) n cos ( 2 n arcsin ( x ) = ( − 1 ) n n ∑ k = 0 n ( − 1 ) k ( n + k − 1 ) ! ( n − k ) ! ( 2 k ) ! 4 k x k , p n ( x ) = ∑ k = 0 n ( − 1 ) n + k n ( n + k − 1 ) ! ( n − k ) ! ( 2 k ) ! 4 k ︸ = a n , k x q k , we get we following series for Dirichlet L-functions, (7) L ( s , χ ) = − 1 1 − χ ( 2 ) 2 1 − s ∑ l = 1 q χ ( l ) ( − 1 ) l ∑ j = 0 n − 1 ( − 1 ) j c n , j ( l + q j ) s + ξ n ( s ) , where c n , j = 1 − ∑ k = 0 j ( − 1 ) k a n , k p n ( − 1 ) = 1 − ∑ k = 0 j α n , k ∑ k = 0 n α n , k , and α n , k = ( n + k − 1 ) ! ( n − k ) ! ( 2 k ) ! 4 k . Note that Borwein (2000) recommended similar approach for the calculation of the Riemann zeta function to make the error as small as possible. Indeed, Chebyshev polynomials T n ( x ), shifted to [ 0 , 1 ] and appropriately normalized maximization of the value p n ( − 1 ) over all polynomials of comparable supremum norm on [ 0 , 1 ].

Let Ω n be a integral random variable associated with Chebyshev polynomials, P ( Ω n = k ) = α n , k ∑ j = 0 n α n , j , 0 ⩽ k ⩽ n . Suppose that F n ( x ) is the cumulative distribution function of the random variable Ω n . It has been proved (cf. Theorem 3.1 in Belovas and Sakalauskas, 2018), that the cumulative distribution function is F n ( x ) = Φ ( x − μ n σ n ) + O ( 1 n ) . By Lemma 2.2 of Belovas and Sakalauskas (2018), the mean and variance are given by μ n = E ( Ω n ) = n 2 + 33 144 2 1 n + O ( 1 n 2 ) , ( σ n ) 2 = V a r ( Ω n ) = 2 8 n + O ( 1 n 2 ) . We now can formulate a result establishing a normal approximation for the coefficients of the series associated with Dirichlet L-functions. Proposition 3.

The following the Gaussian approximation for c n , j coefficients holds: (8) c n , j = 1 − Φ ( j + 1 / 2 − μ n σ n ) + O ( 1 n ) . Note that 1 / 2 is the continuity correction, and σ n = n 2 2 4 + O ( 1 n 5 / 2 ) .

Modern GPUs perform arithmetic instructions much faster when using single-precision (float) instead of double-precision (double) type variables NVIDIA (2023). If the precision using float type variables is good enough (i.e. there is no visible difference between images generated using float and double type variables), this precision for performance trade-off is preferred. However, the EMB-algorithm requires to compute a large sum of values very close to zero. This summation will compound errors that occur when using single-precision variables and produce a result that is not good enough for visualization.

In this Subsection, we aim to determine the domain within which the EMB-algorithm struggles to produce acceptable results. The experiments have been performed on GeForce GTX 1060 (Mobile) graphics processor with 6 GB VRAM and Intel^® Core^TM i7-8750H 2.2 GHz (boosted to 4.0 GHz) processor with 16 GB DDR4 RAM.

The first experiment showed that an algorithm using float type variables is about 22 times faster than an equivalent algorithm using double type variables. However, the accuracy of the computations has dropped dramatically.

The second and third experiments aimed to evaluate the accuracy of the float algorithm. In the second experiment, we showed that the accuracy decreases as t increases. This makes sense because, as t increases, so does n. This means that the sum that gives us the Riemann zeta value contains more operations and, subsequently, more rounding errors. We have also noticed that the highest errors were obtained near the special points. In the third experiment, we further investigated the areas around special points. From the results, we can conclude that within a radius r = 1 / 2 starting from the 25th special point, we can only expect 3 decimal digits of accuracy.

The low accuracy of the float algorithm makes it less suitable for any high-precision Riemann zeta value computations. However, it might still be suitable for visualizations where low precision might be enough. This phenomenon will be investigated further in the next section.

A classic method of visualizing 3D objects is using a ray casting graphics engine. This requires creating a polygon mesh representing an approximation of the object’s shape. A polygon mesh is a set of vertices connected by edges where three edges create a face. In the case of a Mandelbulb, the polygon mesh would have to represent the surface area between the regions that do and do not belong to the Mandelbulb. For this task, we can take inspiration from the 2D Mandelbrot fractal. In Fig. 5, we can see a closeup of the Mandelbrot set where the points that belong to the set are marked with a plus sign and the points that do not are marked with a minus sign. If we were a 2D observer, when placed in the minus region of the image, we could only see the edge of the fractal (the red cells). These red cells can be easily identified if they belong to the Mandelbrot set and have at least one adjacent cell that does not belong to the Mandelbrot set. For a 3D equivalent, we can divide the volume into cubes. Then, we can calculate whether the centre point of the cube belongs to the Mandelbulb and, thus, whether the cube should be visible. This approach would create a blocky polygon mesh. Surface smoothing algorithms could be used to smooth the surface of the polygon. This polygon might be suitable with a high enough resolution and looking from far enough away. However, we would be limited by the mesh’s size if we were to explore the fractal details of a Mandelbulb. The created model is also very rigid, meaning that if we would like to view a Mandelbulb of a different power, we would have to create an entirely new mesh. For these reasons, the ray casting approach is unsuitable for our task. However, it could still be useful for printing a Mandelbulb.

A technique called ray marching is a way to visualize three-dimensional objects without using a polygon mesh of the object. This method instead uses signed distance functions (SDFs) to estimate the distance to an object. An SDF of an object takes a point’s coordinates as an argument and returns the distance to the closest point of the surface of that object. Then we modify a simple ray casting engine to use these SDFs to render an image of an object. SDFs can also be used to calculate surface normals of an object. This means that objects rendered using the ray marching method can be shaded. We can use ray marching to render Mandelbulbs. For that purpose, we need an SDF for a Mandelbulb, which luckily exists. Algorithm 1 (Syntopia, 2011) gives us the SDF for a Mandelbulb.

Zetabulb visualization requires associating the Mandelbulb signed distance function with the Riemann zeta function. When combining a fractal with another object (e.g. shape, fractal, surface), we transform the fractal’s space by the object’s properties and then visualize the fractal in this new space. Thereby, we will combine the Mandelbulb with the surface of the Riemann zeta function.

The procedure is nontrivial. The Mandelbulb has a fractal surface; however, its volume is finite, whereas the Riemann zeta function’s domain is an infinite complex plane. Thus, the question arises: how do we combine a finite volume and an infinite plane? We will solve this problem by adding two more degrees of freedom: movements in the complex plane of the Riemann zeta function. This would allow the Mandelbulb to inherit the properties of a small section of the Riemann zeta functions surface rather than the entire plane. We can implement the new degrees of freedom by adding a complex point c. This point will associate the Mandelbulb with the Riemann zeta function.

The association of the complex plane and the Cartesian coordinate system can be accomplished in two steps. First, we compute the Riemann zeta function value z for an argument related to the point c. Next, we construct a new point η and compute the SDF of a Mandelbulb for this point. An example of this association is presented in Algorithm 2. Figure 6 displays Zetabulbs obtained by several association methods. We will use the same principle when associating Mandelbulbs with Dirichlet L-functions.

We have discovered that algorithms using float type variables perform ∼ 22 times faster than equivalent ones employing double type variables. However, the second experiment showed that this gain comes with a significant loss in accuracy. In this section, we will determine how this loss of accuracy impacts Zetabulb image quality.

Algorithm 3 defines the process of how each pixel is coloured. In total, the Riemann zeta function is called up to maximumSteps times and an additional 6 times when computing a surface normal if the ray hits the Zetabulb. Errors from having non-accurate Riemann zeta values can quickly accumulate and produce incorrect pixel colour values. We will perform an experiment to test how inaccurate the computed values are. First, we will generate ten Zetabulb images with 1 ⩽ d ⩽ 10 decimal digits of accuracy with the double algorithm and one with the float algorithm. We will call the image with 10 decimal digits of accuracy our benchmark and the others our test images. Then, we will calculate the colour accuracies of each image when compared to the benchmarks’ corresponding pixels.

To compute each pixel’s colour accuracies Δ E, we will use the CIEDE2000 algorithm (see Sharma et al. (2005)), an industry-standard, to evaluate how similar colours are, Δ E ∈ [ 0 ; 100 ]. Empirically, these values can be categorized into the following groups:

In our experiment, we compute the percentages of pixels of our test images that fall in each group. By comparing the percentages of the float algorithm to the other test images, we should be able to evaluate the adequate visualization accuracy of the float algorithm.

Earlier, we have established that our implementation of the MB-algorithm will compute the coefficients { ψ n , k } for a maximum imaginary part t max and d decimal digits of accuracy. The same coefficients { ψ n , k } will be used to compute Riemann zeta function values for t < t max . Naturally, these values will be computed with higher decimal accuracy than d. This property of our algorithm makes it problematic to compute the Riemann zeta function values for arguments with different imaginary parts while maintaining exactly the same accuracy. For this reason, as a benchmark, we use values computed using the double algorithm with d = 10 decimal digits of accuracy. Next, by rounding down the benchmark values, we obtain Riemann zeta function values with lower accuracies.

For this experiment, we evaluate the visualization accuracy in two scenarios, with ℑ c = 5 and with ℑ c = 500. We also produce a grayscale image of the float image with highlighted pixels that fall into the groups Δ E ∈ ( 10 ; 50 ] and Δ E ∈ ( 50 ; 100 ]. The computed Δ E values are presented in Tables 3 and 4. Figure 7 presents the generated images for each scenario.

According to Table 3, in Scenario 1 ( ℑ c = 5 ) the effective visualization quality could be achieved between d = 6 and d = 7. However, we can notice an unusually high number of very different pixels (groups Δ E ∈ ( 10 ; 50 ] and Δ E ∈ ( 50 ; 100 ]) for these quality ranges. In Fig. 7(f), we can see these pixels highlighted by the colours red and green, respectively. We can see a small green region in the bottom right of the image. The float algorithm appears to have generated a segment that does not exist in the benchmark image. After further investigation, we could see that the region exists, and the double algorithm could not generate it. Scenario 2 ( ℑ c = 500 ) did not have such anomalies, and its effective precision could be evaluated as between d = 4 and d = 5 in all Δ E categories. As we can see, by increasing ℑ c, we experienced a degradation in image visualization quality.

With the third experiment, we have discovered that the float algorithm is inaccurate when computing the Riemann zeta function values in the neighbourhoods of the special points s k . Scenario 3 aims to investigate Zetabulbs generated near one of these points. Figure 8(a) was generated using the double algorithm with d = 10 decimal digits of accuracy, whereas Fig. 8(b) was generated using the float algorithm. We can see that the float algorithm generated different although similar structures.

In this subsection, we evaluate the real-time Zetabulb visualization performance. We measure the average FPS (frames per second) that our Zetabulb visualization instrument can achieve with float and double algorithms and different GPUs in four experiments, each containing four scenarios lasting for 30 seconds:

The visualization performance results are presented in Table 5. The human eye needs at least 24 FPS to perceive a video as smooth, so this is the value that we will consider satisfactory (coloured blue in Table 5). As discussed before, using double type variables in CUDA computations significantly increases performance time. We can see this effect again when comparing the results of T 1 to T 3 and T 2 to T 4 where the double algorithm could not achieve even a single frame per second.

Scenarios A and B are used as benchmarks. As discussed before, this work required the creation of a custom ray marching graphics engine to test the EMB-algorithm efficiently. In order to show that this engine does not hinder the results, regular Mandelbulb images were generated. Scenarios A and C produce static scenes, whereas B and D produce dynamic scenes. We can see no difference in performance between static and dynamic scene results. This can be viewed as both an advantage and a disadvantage. Our tool does not have a feature that would cache pixel values if no scene parameters change between frames. However, this also means that changing scene parameters between frames has no impact on the performance, so we can create videos of transforming fractals at high frame rates.

The real-time Zetabulb visualization performance results were poor when using the GTX 1060; however, the RTX 3080 produced satisfactory results in T 1 . It should be noted that the RTX 3080 GPU is now over three years old, and faster GPUs exist on the market. The usage of newer GPUs could produce even better performance results. We conclude that the proposed algorithm produces satisfactory results with low t max values and higher-end GPUs.

In this subsection, we will provide Zetabulb images with c values near interesting features of the Riemann zeta function surface. All of the Figs. 10(a) through 10(b) depicting Zetabulbs were generated using Algorithm 3 where the SDF is set to Algorithm 2 using both the Riemann zeta function (which we will call Riemann Zetabulbs) and the Dirichlet L-function (Dirichlet Zetabulbs) with the character χ 5 , 4 (exact value shown in Table 6).

The argument value c for Algorithm 2 has to be chosen carefully. Let us consider ε-neighbourhood of the point c. The algorithm will produce normal Mandelbulbs if | ζ ( c + ε ) | ≈ 0. If | ζ ( c + ε ) has very large values, the algorithm will not produce any fractal structures. Ideally, we want to position the point c where | ζ ( c + ε ) | ∈ ( 0 ; 5 ] and where there is a lot of variance in the nearby values. The Riemann zeta function values can be divided into five regions: the pole, the critical strip, the plateau to the right of the critical strip where ζ ( c ) ≈ 1, the wall to the left of the critical strip where ζ ( c ) values are very large, and the small triangular plateau to the left of the pole where ζ ( c ) ≈ 0. Figure 9 shows all five of these regions. The points that fit the criteria for generating fractal structures lie on the boundaries of all regions and the entire critical strip.

We selected three points to showcase in this work. As the first point, we chose the pole of the Riemann zeta function c = 1. Here, the large values of the pole region completely remove the middle of the Mandelbulb, leaving a bow-shaped structure with holes, which can be seen in Figs. 10(a) and 10(b). For comparison, the Dirichlet L-function does not have a pole, thus the Zetabulb in Fig. 10(c) appears similar to a regular Mandelbulb.

We chose the seventh trivial zero of the Riemann zeta function c = ρ 7 for the second point of interest. The corresponding fractal is visualized in Fig. 11(a). The Riemann zeta function has trivial zeroes at points ρ k = − 2 k, k ∈ N. The trivial zeroes ρ k , with k ∈ [ 1 ; 6 ] lie in the triangular plateau. An example of one such fractal is presented in Fig. 11(b). Whereas, the trivial zeroes with k ∈ [ 8 ; ∞ ) belong to the wall region. This leaves only the point ρ 7 to produce distinct structures when compared to the original Mandelbulb (see Fig. 11) which is calculated using Algorithm 1.

As the final point of interest, we have selected the fourth nontrivial zero of the Riemann zeta function c = γ 4 = 1 / 2 + 30.42487 i. At this point, most of the fractal structure forms a round blob, whereas the top of the fractal is slightly extended, creating a bottle-like shape that can be seen in Figs. 12(a) and 12(b).

•

Two versions of the optimal approach (one using float type variables and one using double type variables) have been compared. The conclusion is that the float algorithm is much faster ( ∼ 22 times) while being less accurate and unsuitable for precise Riemann zeta function computation.

The presented 3D fractal (Zetabulb) images generating computer graphics software tool could be used as a base for SFX and CGI in film-making, as well as the investigation of underlying properties of the Riemann zeta function.