1 Introduction
Since their introduction, Mandelbrot sets have been studied extensively and have captivated both mathematicians and fractal enthusiasts. It was therefore natural to seek a three-dimensional analogue of the Mandelbrot set. A direct three-dimensional counterpart, however, cannot exist in a canonical sense, as there is no true three-dimensional analogue of the complex plane. To overcome this limitation, researchers developed alternative structures based on spherical coordinates, enabling the visualization of a three-dimensional representation of Mandelbrot-type dynamics. These structures became known as Mandelbulbs (cf. Ugwu et al., 2024; Kumar and Bhaskar, 2025).
Beyond their mathematical appeal, Mandelbulbs have demonstrated practical relevance. They have been employed in the creation of digital visual effects for high-budget films, illustrating that these structures are not merely objects of mathematical aesthetics, but also tools with concrete applications in computational art and media.
Recently, we have further developed (see Belovas et al., 2022) an efficient algorithm originally proposed by Borwein (2000) for the computation of the Riemann zeta function. This algorithm is particularly well suited for the evaluation of multiple zeta values (cf. Algorithm 1 in Belovas et al. (2022)). The Mandelbrot set associated with the Riemann zeta function exhibits intricate and highly nontrivial structures in the neighbourhoods of its pole and zeros (see Belovas et al., 2022). These observations naturally lead to the problem of a three-dimensional analogues, namely, Mandelbulbs associated with zeta functions.
How do the resulting fractal structures behave in neighbourhoods of poles and nontrivial zeros of zeta functions? Do such structures possess mathematical and practical significance? The present study aims to advance our understanding of these questions. Specifically, we seek to construct and investigate Mandelbulbs associated with the Riemann zeta function and Dirichlet L-functions, which we refer to as Zetabulbs. To achieve this objective, the following tasks must be undertaken:
-
• implementation of GPU-based (Graphics Processing Unit) versions of efficient algorithms for the computation of zeta functions (modified Borwein’s and Šleževičienė’s algorithms);
-
• development of efficient algorithms for the real-time generation and rendering of Mandelbulb images;
-
• construction of Mandelbulb dynamics derived from zeta functions;
-
• visualization of Zetabulbs and exploration of newly received fractal structures.
2 Calculation of the Riemann Zeta Function
The Riemann zeta-function is a function of a complex variable $s=\sigma +it$ defined by the ordinary Dirichlet series or the Euler product formula,
\[ \zeta (s)={\sum \limits_{n=1}^{\infty }}\frac{1}{{n^{s}}}=\prod \limits_{p\hspace{0.2778em}\text{prime}}{\bigg(1-\frac{1}{{p^{s}}}\bigg)^{-1}},\hspace{1em}\mathrm{\Re }s\gt 1,\]
and by analytic continuation elsewhere except for a simple pole at $s=1$. It satisfies the functional equation
implying that $\zeta (s)$ has simple zeros at $s=-2n$, $n\in \mathbb{N}$, known as the trivial zeros. The Riemann hypothesis states that all other zeroes, called non-trivial, lie on the critical line $\mathrm{\Re }s=1/2$.To investigate the Riemann zeta function numerically, one needs an algorithm enabling him to effectively compute multiple values over the complex plane. Moreover, the algorithm must be well-suited for parallel CPU (Central Processing Unit) and GPU computations. We have employed modifications of Borwein’s efficient algorithm. The idea of the method comes from the alternating series convergence and properties of Chebyshev polynomials (see Borwein, 2000). It is important to note that the algorithm is nearly optimal in the sense that there is no sequence of n-term exponential polynomials that can converge to $\zeta (s)$ very much faster than those of the algorithm (cf. Theorem 3.1 in Borwein, 2000).
2.1 MB-Algorithm
The modified Borwein’s algorithm (MB-algorithm) applies to complex numbers with $\mathrm{\Re }s\geqslant 1/2$ and arbitrary $\mathrm{\Im }s$. In this modification (see Belovas et al., 2022), the Riemann zeta function is represented by the alternating series
Here
and
(1)
\[ \zeta (s)=\frac{1}{1-{2^{1-s}}}{\sum \limits_{k=0}^{n-1}}\frac{{(-1)^{k}}{\psi _{n,k}}}{{(k+1)^{s}}}+{R_{n}}(s).\]
\[ {l_{\max }}=\text{arg}\underset{0\leqslant k\leqslant n}{\max }\frac{(n+k-1)!{4^{k}}}{(n-k)!(2k)!},\hspace{1em}n\in \mathbb{N},\hspace{2.5pt}0\leqslant k\leqslant n,\]
while
\[\begin{aligned}{}& {H_{l}}={H_{l-1}}+\exp \big({T_{l}}-{T_{{l_{\max }}}}+(l-{l_{\max }})\log 4\big),\\ {} & {H_{0}}=\exp ({T_{0}}-{T_{{l_{\max }}}}-{l_{\max }}\log 4),\\ {} & {T_{l}}={T_{l-1}}+\log \frac{(n-l+1)(n+l-1)}{(2l-1)(2l)},\\ {} & {T_{0}}=-\log n,\hspace{1em}1\leqslant l\leqslant n.\end{aligned}\]
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
here $|s-{s_{k}}|\geqslant \varepsilon \gt 0$ and ${s_{k}}=1+2i\pi k/\log 2$.
(2)
\[ n=\bigg\lceil \frac{\pi t}{2\log (3+\sqrt{8})}+\frac{\log 10}{\log (3+\sqrt{8})}d+\frac{\log \frac{2}{1-{2^{-\varepsilon }}}}{\log (3+\sqrt{8})}\bigg\rceil ,\]2.2 NA-Algorithm
The limit theorem for the coefficients of MB-series enables us to derive the following normal approximation for the coefficients ${\psi _{n,k}}$ of the series (1) (cf. Belovas et al., 2022)
Note that for fixed $\mathrm{\Re }s$ and d, we have ${k_{0}}\sim 0.63\mathrm{\Im }s$.
\[ {\psi _{n,k}}=1-\Phi \bigg(\frac{k-{\mu _{n}}}{{\sigma _{n}}}\bigg)+O\bigg(\frac{1}{\sqrt{n}}\bigg),\hspace{1em}{\mu _{n}}=\frac{n}{\sqrt{2}},\hspace{1em}{\sigma _{n}}=\frac{\sqrt{n}}{\sqrt[4]{32}}.\]
Here $\Phi (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}}=\lceil {\mu _{n}}+z{\sigma _{n}}\rceil ,\hspace{1em}z={\Phi ^{-1}}\big(1-{10^{-d}}\big),\hspace{2.5pt}n\gt {10^{3}}.\]2.3 EMB-Algorithm
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
The empirical formula (4) has been verified for $0\leqslant t\leqslant {10^{4}}$.
(4)
\[ \begin{aligned}{}& n=\lceil at+b\sqrt{t}+c\rceil ,\hspace{1em}a=0.451,\\ {} & b=1.407\sqrt{d}-0.245,\hspace{1em}c=0.371d+0.195.\end{aligned}\]3 Calculation of the Dirichlet L-Functions
Let $s\in \mathbb{C}$ and $\mathrm{\Re }s\gt 1$, then Dirichlet L-series are defined as follows
\[ L(s,\chi )={\sum \limits_{n=1}^{\infty }}\frac{\chi (n)}{{n^{s}}},\hspace{1em}\mathrm{\Re }s\gt 1.\]
Here χ is a Dirichlet character modulo q. By analytic continuation, this function extends to a meromorphic function on the whole complex plane and is known as a Dirichlet L-function. The Dirichlet L-function is a generalization of the Riemann zeta function and thus has similar properties. If χ is principal, the corresponding Dirichlet L-function has a simple pole at $s=1$; otherwise, it is holomorphic on the whole complex plane. The Generalized Riemann hypothesis asserts that all nontrivial zeros of Dirichlet L-function lie on the critical line $\mathrm{\Re }s=1/2$.Šleževičiene (2004), building upon Borwein’s technique, proposed several algorithms (MS-algorithms) for efficient computing of Dirichlet L-functions.The following two results (Propositions 1 and 2), established in Šleževičiene (2004) (see Theorems 4 and 5, respectively), yield algorithms for computing Dirichlet L-functions attached to a character with odd modulus.
Proposition 1.
Let q be odd, χ mod q be a non-principal character and let
be an arbitrary polynomial of degree $nq$ which does not vanish at $-1$. Further, let
Then for $\sigma \gt 0$,
where
Proposition 2.
Let q be odd. For $\sigma \gt 0$,
where
(6)
\[ L(s,\chi )=\frac{1}{{2^{n}}(1-\chi (2){2^{1-s}})}{\sum \limits_{a=1}^{q}}{(-1)^{a-1}}\chi (a){\sum \limits_{j=0}^{2n-1}}\frac{{e_{j}}}{{(a+qj)^{s}}}+{\gamma _{n}}(s),\]
\[ {e_{j}}={(-1)^{j}}\Bigg({\sum \limits_{k=0}^{j-n}}\left(\genfrac{}{}{0pt}{}{n}{k}\right)-{2^{n}}\Bigg).\]
If $\varphi (m)$ stands for the Euler totient function, then
In the next section, we introduce a modification of Šleževičienė’s efficient algorithm for Dirichlet L-functions that applies a normal approximation to the coefficients of the corresponding series.
3.1 NA-Modification of MS-Algorithm for Dirichlet L-Functions
Associating polynomials ${p_{n}}(x)$ in Proposition 1 with Chebyshev polynomials
where
\[\begin{aligned}{}{T_{n}}(x)& ={(-1)^{n}}\cos (2n\arcsin (\sqrt{x})\\ {} & ={(-1)^{n}}n{\sum \limits_{k=0}^{n}}{(-1)^{k}}\frac{(n+k-1)!}{(n-k)!(2k)!}{4^{k}}{x^{k}},\\ {} {p_{n}}(x)& ={\sum \limits_{k=0}^{n}}\underset{={a_{n,k}}}{\underbrace{{(-1)^{n+k}}n\frac{(n+k-1)!}{(n-k)!(2k)!}{4^{k}}}}{x^{qk}},\end{aligned}\]
we get we following series for Dirichlet L-functions,
(7)
\[ L(s,\chi )=\frac{-1}{1-\chi (2){2^{1-s}}}{\sum \limits_{l=1}^{q}}\chi (l){(-1)^{l}}{\sum \limits_{j=0}^{n-1}}\frac{{(-1)^{j}}{c_{n,j}}}{{(l+qj)^{s}}}+{\xi _{n}}(s),\]
\[ {c_{n,j}}=1-\frac{{\textstyle\textstyle\sum _{k=0}^{j}}{(-1)^{k}}{a_{n,k}}}{{p_{n}}(-1)}=1-\frac{{\textstyle\textstyle\sum _{k=0}^{j}}{\alpha _{n,k}}}{{\textstyle\textstyle\sum _{k=0}^{n}}{\alpha _{n,k}}},\]
and
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 ${\Omega _{n}}$ be a integral random variable associated with Chebyshev polynomials,
\[ P({\Omega _{n}}=k)=\frac{{\alpha _{n,k}}}{{\textstyle\textstyle\sum _{j=0}^{n}}{\alpha _{n,j}}},\hspace{1em}0\leqslant k\leqslant n.\]
Suppose that ${F_{n}}(x)$ is the cumulative distribution function of the random variable ${\Omega _{n}}$. It has been proved (cf. Theorem 3.1 in Belovas and Sakalauskas, 2018), that the cumulative distribution function is
\[ {F_{n}}(x)=\Phi \bigg(\frac{x-{\mu _{n}}}{{\sigma _{n}}}\bigg)+O\bigg(\frac{1}{\sqrt{n}}\bigg).\]
By Lemma 2.2 of Belovas and Sakalauskas (2018), the mean and variance are given by
\[\begin{aligned}{}& {\mu _{n}}=E({\Omega _{n}})=\frac{n}{\sqrt{2}}+\frac{33}{144\sqrt{2}}\frac{1}{n}+O\bigg(\frac{1}{{n^{2}}}\bigg),\\ {} & {({\sigma _{n}})^{2}}=Var({\Omega _{n}})=\frac{\sqrt{2}}{8}n+O\bigg(\frac{1}{{n^{2}}}\bigg).\end{aligned}\]
We now can formulate a result establishing a normal approximation for the coefficients of the series associated with Dirichlet L-functions. Proposition 3.
Proof.
The proposition is an immediate consequence of the proof of Theorem 2 in Belovas et al. (2022). Besides,
□
4 Parallel Algorithms for Calculation of the Riemann Zeta Function
Two approaches to parallelize the MB-algorithm can be used: optimize for specific or multiple value computation speed. Each pixel of an image of a Zetabulb will require multiple computations of the Riemann zeta function. We aim to generate at least 24 frames of HD: 720p images (a total of $720\times 1280=921600$ pixels) per second. For this purpose, optimizing for multiple-value computation speed is more beneficial.
Personal computers usually have two processing units: a central processing unit (CPU) and a graphics processing unit (GPU). Both processing units can have multiple cores; however, they are specialized for different tasks. CPUs usually have a low amount (up to 24) of high-speed (up to 6GHz) cores and are used for tasks that are hard to parallelize. Whereas, GPUs usually have a high amount (up to 16000) of slower (up to 2.5GHz) cores and are used for tasks that are easy to parallelize. CUDA NVIDIA (2026) is a parallel computing platform that allows one to perform computations on both the CPU and the GPU and manages the data transfer operations between RAM (Random Access Memory) and VRAM (Video Random Access Memory). This platform is the optimal approach for real-time image generation.
The MB-method can be separated into two logical parts: initialization and computation. During the first part, we compute n coefficients ψ required to compute the Riemann zeta function with the maximum imaginary part ${t_{\max }}$ and d decimal digits of accuracy. During the second part, we use these coefficients to compute the actual value of the Riemann zeta function. Note that although both parts are linear, the coefficients computed during the initialization part can be reused when computing Riemann zeta values for arguments $t\lt {t_{\max }}$. This property of the MB-method makes it suitable for computing multiple values at once.
A data transfer diagram for a CUDA implementation of the MB-method can be seen in Fig. 1. Here, the data flow of the initialization part is indicated by blue arrows and the computation part by red arrows. During the initialization part, the coefficients ψ are computed using a CPU. After that, they are transferred from RAM to VRAM for the remainder of the program’s runtime. Note that it might be more efficient to perform the initialization part on the GPU if the CPU clock speed is low. Then, whenever the CPU requests a Riemann zeta value to be computed, the GPU uses ψ to get the value and transfers it from VRAM back to RAM. An implementation of this algorithm can be found in Kuzma (2025).
Fig. 1
Data transfer diagram for the MB-method. The blue arrows indicate the initialization part’s data transfer flow. The red arrows indicate the computation part’s data transfer flow.
4.1 Numerical Experiments
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® CoreTM i7-8750H 2.2 GHz (boosted to 4.0 GHz) processor with 16 GB DDR4 RAM.
4.1.1 First Experiment
The first experiment aims to compare the speed and accuracy of the float and double algorithms. For this experiment, we divide the complex plane into ${n_{f}}=1000$ matrices (frames) consisting of $1280\times 720$ uniformly located points
where ${k_{p}}={10^{2(p-1)}}$ is the minimum imaginary value of a frame, ${h_{f}}=\frac{{k_{p+1}}-{k_{p}}}{{n_{f}}}$ is the frame height and $0\lt j\leqslant {n_{f}}$. We compute the Riemann zeta function values once with an algorithm that uses double precision variables and once with an algorithm that uses float precision variables. We find the average mean square error and the maximum mean square error of the computed values. Using both methods, we also compute the average times to compute a frame. The results of this experiment are shown in Table 1.
(9)
\[ {F_{p,j}}=\underset{\sigma }{\underbrace{(0.5,1.5)}}\times \underset{t}{\underbrace{\big({k_{p}}+j{h_{f}},{k_{p}}+(j+1){h_{f}}\big)}},\]Table 1
First experiment results. $AE$ stands for average root mean square error between float and double frames, $ME$ stands for maximum square error between the float and double frames, ${\tau _{\textit{double}}}$ – average time to compute a frame using the double algorithm in milliseconds, ${\tau _{\textit{float}}}$ – average time to compute a frame using the float algorithm.
| p | $AE$ | $ME$ | ${\tau _{\textit{double}}}$ | ${\tau _{\textit{float}}}$ |
| 2 | $2.43\cdot {10^{-5}}$ | $9.25\cdot {10^{-3}}$ | 201.4 | 9.7 |
| 3 | $9.38\cdot {10^{-4}}$ | $8.35\cdot {10^{-2}}$ | 32729.3 | 1459.4 |
4.1.2 Second Experiment
In the second experiment, we aim to evaluate the float algorithm’s precision on the critical line. We aim to find the maximum value ${t_{\min }}$ with which the float algorithm no longer produces d decimal digits of accuracy. In Belovas et al. (2022), we have shown that the most computationally expensive region of the complex plane (excluding the special points ${s_{k}}=1+2\pi ki/\log 2$) for the MB-algorithm is the critical line.
We compute the Riemann zeta function values using both the float and the double algorithms and compute the error for each of the obtained values $E=|\zeta {({F_{k,j}})_{\textit{float}}}-\zeta {({F_{k,j}})_{\textit{double}}}|$. Then, we find the first value t with which the float algorithm no longer provides d decimal digits of accuracy. Here
The experiment results are presented in Fig. 2. As we can see, the error E increases with t. In Table 2, we have listed the largest imaginary part of the complex number $1/2+it$ that produces d decimal digits of accuracy. For example, we can expect that every t value lower than 72.52 will be computed with at least 3 decimal digits of accuracy using the float algorithm.
4.1.3 Third Experiment
As already noted, the special points ${s_{k}}=1+2\pi ki/\log 2$ are the most challenging to compute using the MB-algorithm. In the third experiment, we sought to estimate the radius r of a circle within which we could expect d decimal digits of accuracy could be expected for the first 200 special points ${s_{k}}$. To this end, we constructed a ${M_{2000\times 2000}}$ uniform square grid of complex points centred at each special point with $\sigma \in [1/2;3/2]$ and $t\in [2\pi k/\log 2-1/2;2\pi k/\log 2+1/2]$. We repeated the experiment using ${2^{22}}$ points of the Sobol sequence. The comparison showed that the resulting error distributions were qualitatively similar to those obtained with the uniform grid. The corresponding results are presented in Fig. 4. Next, we calculated the Riemann zeta function for each of these points using the double (with $d=10$ decimal digits of accuracy) and float algorithms and created a new $2000\times 2000$ grid of log-errors
\[ {\hat{d}_{n,j}}=\big\lfloor -\lg \big|\zeta {({M_{n,j}})_{\textit{float}}}-\zeta {({M_{n,j}})_{\textit{double}}}\big|\big\rfloor ,\hspace{1em}1\leqslant n,j\leqslant 2000.\]
Note that ${\hat{d}_{n,j}}\in \mathbb{N}$ and ${\hat{d}_{n,j}}\in [1;10]$; however, we have obtained $\max {\hat{d}_{n,j}}=4$. A visual example of this grid for the 1st, 10th, 100th and 1000th special points can be seen in Fig. 3.Next, we find the maximum distance for each fixed ${\hat{d}_{n,j}}=l$,
\[ {r_{k,l}}=\underset{{\hat{d}_{n,j}}=l}{\max }|{M_{n,j}}-{s_{k}}|,\hspace{1em}1\leqslant l\leqslant 4,\hspace{1.5pt}1\leqslant k\leqslant 200.\]
The experiment results are presented in Fig. 4. The maximal ${r_{k,l}}$ value, ${r_{k,l}}=1/\sqrt{2}$, is reached at the corner of the square grid (see Fig. 4). We can see that as a special point’s t value increases, so does the area around which Riemann zeta values are computed with l decimal digits of accuracy. It should also be noted that this precision area is not perfectly circular or uniform. This means we can still expect to compute higher precision Riemann zeta values within these areas.4.2 Parallelization Methods Comparison
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/\sqrt{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.
5 Visualization of 3D Fractals
Before rendering a Mandelbulb, recall how a two-dimensional Mandelbrot set is generated. The set consists of complex numbers c for which the function ${f_{c}}(z)={z^{2}}+c$ does not diverge when iterated from $z=0$. The pixels representing these points are marked with a plus sign in Fig. 5. The other points are usually coloured according to the iteration number, at which the value exceeds some preset threshold. For a 2D fractal, this approach works effectively: pixels can be coloured on a plane, and the observer perceives the overall fractal structure from above (cf. Belovas et al., 2022). In contrast, while considering a three-dimensional analogue, the observer is placed within the diverging or non-diverging regions. In this context, the observer primarily perceives the boundary between the two regions rather than the regions themselves. Consequently, the traditional method of colouring pixels based on iteration count is insufficient for 3D visualization.
In this section, we outline two approaches to the visualization of Mandelbulbs, addressing these challenges and enabling the rendering of three-dimensional fractal structures.
Fig. 5
Closeup segment of the Mandelbrot set. The cells marked with the plus sign belong to the Mandelbrot set, and those with the minus sign do not. The striped cells represent the edge of the Mandelbrot fractal.
5.1 Ray Casting Approach
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.
5.2 Ray Marching Approach
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.
6 Zetabulb Visualization
As it has been concluded, we use the ray marching technique to render Zetabulbs. Several free ray marching engines that can visualize Mandelbulbs are available. However, none allow efficient CUDA utilization for the Riemann zeta function calculation. Therefore, we will create a basic ray marching engine from scratch. The engine’s source code and the MB-algorithm’s implementation can be found in Kuzma (2025). The new instrument will exhibit the following features:
6.1 Zetabulb Signed Distance Function
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.
Algorithm 2
Signed distance function of a Zetabulb. Here η – the three-dimensional coordinates of the origin point, p – the power of the Mandelbulb, k – the maximum number of iterations, c – the complex number used to associate the Mandelbulb with the real and imaginary surfaces of the Riemann zeta function.
6.2 Visualization Quality
We have discovered that algorithms using float type variables perform $\sim 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.
Fig. 6
Zetabulbs constructed using different association methods with identical c values. Values z and η correspond to the rows 2 and 3 of Algorithm 2.
Algorithm 3
Ray marching function used for the custom rendering tool. This function fires a ray from $\textit{position}$ in $\textit{direction}$ until it hits or misses the object. Returns a colour for a pixel. It is assumed that the centre of the object defined by the SDF is at coordinates (0, 0, 0).
Algorithm 3 defines the process of how each pixel is coloured. In total, the Riemann zeta function is called up to $\textit{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\leqslant d\leqslant 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 $\Delta E$, we will use the CIEDE2000 algorithm (see Sharma et al. (2005)), an industry-standard, to evaluate how similar colours are, $\Delta E\in [0;100]$. Empirically, these values can be categorized into the following groups:
-
• $\Delta E=0$. The colours are the same.
-
• $\Delta E\in (0;1]$. The colour differences are not perceptible by the human eye.
-
• $\Delta E\in (1;2]$. Perceptible through close observation.
-
• $\Delta E\in (2;10]$. Perceptible at a glance.
-
• $\Delta E\in (10;50]$. Colours are more similar than different.
-
• $\Delta E\in (50;100]$. Colours are more different than similar.
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 $\{{\psi _{n,k}}\}$ for a maximum imaginary part ${t_{\max }}$ and d decimal digits of accuracy. The same coefficients $\{{\psi _{n,k}}\}$ will be used to compute Riemann zeta function values for $t\lt {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.
Table 3
E values of the test images for Scenario 1. d – the decimal digits of accuracy with which an image was generated, and f stands for the float algorithm. Other columns lists the percentage of pixels that fall into each $\triangle E$ category when compared with the benchmark image.
| $\Delta E\in $ | ||||||
| d | 0 | $(0;1]$ | $(1;2]$ | $(2;10]$ | $(10;50]$ | $(50;100]$ |
| 1 | 31.93 | 1.68 | 1.52 | 11.47 | 39.96 | 13.44 |
| 2 | 35.04 | 3.71 | 3.15 | 14.13 | 33.97 | 10.00 |
| 3 | 43.39 | 5.58 | 3.34 | 14.29 | 27.59 | 5.80 |
| 4 | 59.41 | 8.83 | 4.74 | 13.65 | 11.50 | 1.86 |
| 5 | 72.65 | 15.18 | 4.49 | 5.42 | 1.96 | 0.29 |
| 6 | 86.45 | 10.94 | 1.10 | 1.08 | 0.38 | 0.05 |
| 7 | 95.02 | 4.21 | 0.30 | 0.33 | 0.12 | 0.02 |
| 8 | 98.55 | 1.23 | 0.08 | 0.09 | 0.04 | 0.00 |
| 9 | 99.81 | 0.17 | 0.01 | 0.01 | 0.00 | 0.00 |
| f | 90.32 | 7.89 | 0.63 | 0.67 | 0.32 | 0.15 |
Table 4
$\triangle E$ values of the test images for Scenario 2. d – the decimal digits of accuracy with which an image was generated, and f stands for the float algorithm. Other columns lists the percentage of pixels that fall into each $\triangle E$ category when compared with the benchmark image.
| $\Delta E\in $ | ||||||
| d | 0 | $(0;1]$ | $(1;2]$ | $(2;10]$ | $(10;50]$ | $(50;100]$ |
| 1 | 11.01 | 2.18 | 8.17 | 15.46 | 59.55 | 9.53 |
| 2 | 16.99 | 4.39 | 2.26 | 17.41 | 51.43 | 6.36 |
| 3 | 23.12 | 5.88 | 3.42 | 21.24 | 41.39 | 3.95 |
| 4 | 37.27 | 17.01 | 4.42 | 22.20 | 12.14 | 1.04 |
| 5 | 54.97 | 32.58 | 10.33 | 4.95 | 1.89 | 0.19 |
| 6 | 83.61 | 14.33 | 5.42 | 0.83 | 0.37 | 0.05 |
| 7 | 94.82 | 4.57 | 0.82 | 0.23 | 0.13 | 0.02 |
| 8 | 98.42 | 1.38 | 0.22 | 0.07 | 0.04 | 0.01 |
| 9 | 99.79 | 0.18 | 0.07 | 0.01 | 0.01 | 0.00 |
| f | 47.81 | 30.78 | 0.01 | 8.86 | 3.89 | 0.49 |
For this experiment, we evaluate the visualization accuracy in two scenarios, with $\mathrm{\Im }c=5$ and with $\mathrm{\Im }c=500$. We also produce a grayscale image of the float image with highlighted pixels that fall into the groups $\Delta E\in (10;50]$ and $\Delta E\in (50;100]$. The computed $\Delta 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 $(\mathrm{\Im }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 $\Delta E\in (10;50]$ and $\Delta E\in (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 $(\mathrm{\Im }c=500)$ did not have such anomalies, and its effective precision could be evaluated as between $d=4$ and $d=5$ in all $\Delta E$ categories. As we can see, by increasing $\mathrm{\Im }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.
Fig. 7
Zetabulb images with different Riemann zeta value precisions for Scenario 1 $(\mathrm{\Im }c=5)$. Subfigures (a) through (d) were generated with the double algorithm and lower precision. Subfigure (e) was generated with the float algorithm. Subfigure (f) is a grayscale version of (e) with pixels that fall into the category $\Delta E\in (10;50]$ highlighted red and $\Delta E\in (50;100]$ highlighted green.
6.3 Visualization Performance
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:
Table 5
Four tests of visualization performance (in FPS). First – float type variables with ${t_{\max }}=50$, second – float type with ${t_{\max }}=500$, third – double type with ${t_{\max }}=50$, fourth – double type with ${t_{\max }}=500$. The required memory in bytes to store the coefficients in VRAM is given for each of the tests. All the tests were performed on two GPUs: RTX 3080 and GTX 1060 Mobile. Values marked in blue denote satisfactory results.
| Variable | float type | double type | ||||||
| Test | ${T_{1}}$ | ${T_{2}}$ | ${T_{3}}$ | ${T_{4}}$ | ||||
| t | ${t_{\max }}=50$ | ${t_{\max }}=500$ | ${t_{\max }}=50$ | ${t_{\max }}=500$ | ||||
| Memory, bytes | 192 | 384 | 1200 | 2400 | ||||
| Scenario | RTX | GTX | RTX | GTX | RTX | GTX | RTX | GTX |
| A | 129.71 | 21.31 | 126.71 | 21.10 | 7.12 | 1.97 | 7.17 | 1.98 |
| B | 137.89 | 22.03 | 142.14 | 21.30 | 6.68 | 1.94 | 7.44 | 1.96 |
| C | 27.35 | 3.99 | 5.78 | 0.97 | 0.11 | 0.03 | 0.02 | 0.005 |
| D | 27.26 | 3.87 | 5.60 | 0.87 | 0.11 | 0.02 | 0.02 | 0.005 |
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.
6.4 Visualization Results
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 ${\chi _{5,4}}$ (exact value shown in Table 6).
Table 6
Dirichlet character ${\chi _{5,4}}$. Note that ${\chi _{5,4}}(-1)={\chi _{5,4}}(4)=1$.
| n | 1 | 2 | 3 | 4 |
| ${\chi _{5,4}}(n)$ | 1 | −1 | −1 | 1 |
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 $|\zeta (c+\varepsilon )|\approx 0$. If $|\zeta (c+\varepsilon )$ has very large values, the algorithm will not produce any fractal structures. Ideally, we want to position the point c where $|\zeta (c+\varepsilon )|\in (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 $\zeta (c)\approx 1$, the wall to the left of the critical strip where $\zeta (c)$ values are very large, and the small triangular plateau to the left of the pole where $\zeta (c)\approx 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.
Fig. 9
The real (yellow) and the imaginary (red) parts of the Riemann zeta function and the five principal regions: the pole neighbourhood (1), the critical strip (2), the plateau to the right of the critical strip (3), the wall to the left of the critical strip (4), and the triangular plateau (5).
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.
Fig. 10
(a) Riemann Zetabulb centered on the pole. (b) Side view of the same Zetabulb. (c) Dirichlet Zetabulb centered on the pole.
We chose the seventh trivial zero of the Riemann zeta function $c={\rho _{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 ${\rho _{k}}=-2k$, $k\in \mathbb{N}$. The trivial zeroes ${\rho _{k}}$, with $k\in [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\in [8;\infty )$ belong to the wall region. This leaves only the point ${\rho _{7}}$ to produce distinct structures when compared to the original Mandelbulb (see Fig. 11) which is calculated using Algorithm 1.
Fig. 11
(a) Riemann Zetabulb centred on the 7th trivial zero. (b) Riemann Zetabulb centred on the 5th trivial zero. (c) Original Mandelbulb.
As the final point of interest, we have selected the fourth nontrivial zero of the Riemann zeta function $c={\gamma _{4}}=1/2+30.42487i$. 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).
Fig. 12
(a) Riemann Zetabulb centred on the fourth non-trivial zero $c={\gamma _{4}}$. Note a bottle-like shape. (b) Closeup image of the bottleneck of the previous Zetabulb. (c) and (d) Dirichlet Zetabulbs centred on the fourth non-trivial zero of the Riemann zeta function. (e) and (f) Dirichlet Zetabulbs centred on the fourth non-trivial zero of the Dirichlet L-function ($c=1/2+16.03382i$).
7 Conclusions
This work aimed to create efficient algorithms, which would be the base of a software instrument for visualizing Mandelbulb fractals associated with the Riemann zeta function and the Dirichlet L-function (Zetabulbs). To the author’s knowledge, this is the first investigation and application of Zetabulbs. To achieve this goal, parallel technologies have been applied to speed up MB-algorithm and technologies designed for fractal visualization and application have been studied and implemented.
7.1 MB-Algorithm Parallelization Conclusions
-
• 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 ($\sim 22$ times) while being less accurate and unsuitable for precise Riemann zeta function computation.
-
• Next, we have noticed that the float algorithm accuracy correlates with the imaginary part of the argument of the Riemann zeta function. Moreover, the highest errors have been obtained in the ε-neighbourhoods of the special points. Hence, the float algorithm should not be used for any precise Riemann zeta function calculations on these domains.
7.2 Zetabulb Visualizations Conclusions
-
• Six different signed distance functions (SDFs) to produce Zetabulb images have been proposed. All the proposed functions have produced satisfactory results in different but still aesthetically appealing fractal structures.
-
• The problem of sufficient (practical) accuracy of the float algorithm compared to double algorithm for Zetabulb visualization has been considered. The conclusion is that the float algorithm is applicable for Zetabulb visualization (for non-scientific, e.g. cinematographic purposes) despite a certain accuracy deterioration with an increase of the imaginary part of the argument.
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.
