Deriving the Error Term Bound Coefficient in the Prime Number Theorem via a Sieve-Based Algorithm

1. Introduction

A prime number (or just prime) is a natural number greater than 1 that is not a product of two smaller natural numbers. Primes are central in number theory because of the fundamental theorem of arithmetic: every natural number greater than 1 is either a prime itself or can be factorized as a product of primes that is unique up to their order. There are infinitely many primes, as demonstrated by Euclid around 300 BC using a very elegant and simple proof.

But how many prime numbers are there up to a given number? Consider the prime-counting function denoted by $\pi \left(x\right)$, which is defined as the number of prime numbers less than or equal to x for any real number x. In 1849, in a letter to the German astronomer Johann Franz Encke [1], Carl Friedrich Gauss made a note suggesting that the number of prime numbers less than or equal to x is:

$$\pi \left(x\right)\sim {\displaystyle \frac{x}{log\left(x\right)}}$$

(1)

The Formula (1) later became known as the Prime Number Theorem (PNT). Legendre conjectured in 1796 that there exists a constant B such that $\pi \left(x\right)$ satisfies:

$$B=\underset{x\to \infty }{lim}\left(logx-{\displaystyle \frac{x}{\pi \left(x\right)}}\right)$$

(2)

B is known as Legendre’s constant. Regardless of its exact value, the existence of this constant implies the Prime Number Theorem. In 1849 Pafnuty Chebyshev proved that if the limit B exists, it must be equal to 1.

The PNT then states that ${\displaystyle \frac{x}{log\left(x\right)}}$ is a good approximation for $\pi \left(x\right)$ (where log here means the natural logarithm), in the sense that the limit of the quotient of the two functions $\pi \left(x\right)$ and ${\displaystyle \frac{x}{logx}}$, as x increases, is 1:

$$\underset{x\to \infty }{lim}{\displaystyle \frac{\pi \left(x\right)}{\left[\frac{x}{logx}\right]}}=1$$

(3)

In a handwritten note on a reprint of his 1838 paper “Sur l’usage des séries infinies dans la théorie des nombres” [2], which he mailed to Gauss, Gustav Lejeune Dirichlet conjectured (under a slightly different form appealing to a series rather than an integral) that an even better approximation to $\pi \left(x\right)$ is given by the offset logarithmic integral function $\mathrm{Li}\left(x\right)$, or the Eulerian logarithmic integral, defined as:

$$\mathrm{Li}\left(x\right)={\int }_2^x{\displaystyle \frac{dt}{lnt}}=\mathrm{li}\left(x\right)-\mathrm{li}\left(2\right)$$

(4)

where $\mathrm{li}\left(x\right)$ is defined as:

$$\mathrm{li}\left(x\right)={\int }_0^x{\displaystyle \frac{dt}{lnt}}$$

(5)

Based on this, the PNT can also be written as:

$$\pi \left(x\right)\sim \mathrm{Li}\left(x\right)$$

(6)

In 1859, Bernhard Riemann in a short but monumental paper titled “Über die Anzahl der Primzahlen unter einer gegebenen Größe” (“On the Number of Primes Less Than A Given Magnitude”) introduced the Riemann zeta function, defined for complex numbers $s=t+i\sigma$ at $\sigma >1$ as (7):

$$\zeta \left(s\right)=\sum _{n=1}^{\infty }{\displaystyle \frac{1}{n^s}}$$

(7)

This function had already been studied by Euler, but for real values of s, and Riemann extended its domain to complex numbers and explored its analytic properties. Crucially, he showed that the distribution of prime numbers is intimately connected to the behavior of the zeros of the zeta function (the values of s for which $\zeta \left(s\right)=0$). In particular, Euler discovered in 1737 the connection between the zeta function and prime numbers by proving that:

$$\sum _{i=1}^{\infty }{\displaystyle \frac{1}{n^s}}=\prod _p{\displaystyle \frac{1}{1-p^{-s}}}$$

(8)

In 1896, Jacques Hadamard and Charles Jean de la Vallée-Poussin independently proved the PNT using complex analysis and properties of the zeta function, specifically, by showing that all non-trivial zeros of the zeta function lie within the critical strip $(0\le \mathrm{Re}(s)\le 1)$. Riemann himself proposed that all non-trivial zeros of the zeta function lie on the critical line $\mathrm{Re}\left(s\right)={\displaystyle \frac{1}{2}}$. This bold conjecture, now known as the Riemann Hypothesis (RH), remains one of the most important unsolved problems in mathematics. In 1914, Godfrey Harold proved that $\zeta ({\displaystyle \frac{1}{2}}+i\sigma )$ has infinitely many real zeros. More than 10 trillion (i.e. ${10}^{13}$) non-trivial zeros of the Riemann zeta function have been verified, and all of them lie on the critical line [3].

The connection between the Riemann zeta function and the prime-counting function $\pi \left(x\right)$ is one reason the Riemann hypothesis has considerable importance in number theory: if established, it would yield a far better estimate of the error involved in the prime number theorem. More specifically, Helge von Koch showed in 1901 that if the RH is true, the error term in the PNT can be expressed as:

$$\pi \left(x\right)=\mathrm{Li}\left(\mathrm{x}\right)+O(\sqrt{x}logx)$$

(9)

where O is a big O notation. This is in fact equivalent to the RH.

The constant involved in the big O notation was estimated in 1976 by Lowell Schoenfeld [4], assuming the RH is true:

$$|\pi \left(x\right)-\mathrm{li}\left(x\right)|<{\displaystyle \frac{\sqrt{x}logx}{8\pi }}$$

(10)

for all $x\ge 2657$. For the Chebyshev prime-counting function $\psi$, which is determined as:

$$\psi \left(x\right)=\sum _{k\in N}\sum _{p^k\le x}logp$$

(11)

a similar bound to (10) can also be derived:

$$|\psi \left(x\right)-x|<{\displaystyle \frac{\sqrt{x}{(logx)}^2}{8\pi }}$$

(12)

The expression for $\pi \left(x\right)$ (10) is connected with the expression for $\psi \left(x\right)$ (12) by the factor $logx$

From the formulas above, it can be seen how strongly the assumption of the truth of the Riemann Hypothesis influences the Prime Number Theorem, particularly in estimating the error term in the distribution. If RH is true, then the error term in the PNT is bounded by the factor $\sqrt{x}logx$.

In this article, an approach to decomposition of $\pi \left(N\right)$ into a dominant term and an error term based on sequential elimination of the composite numbers up to number N is described. This formulation offers a new lens for interpreting the prime distribution and invites further exploration into the analytic and statistical nature of these components. Connections to classical approximations, such as Li(x) and the corresponding error term, and potential links to the zeros of the Riemann zeta function, are also briefly discussed.

2. Estimation of the Prime-Counting Function by Applying of Composite Numbers Elimination

Start by considering the Sieve of Eratosthenes algorithm. It is a simple algorithm to find all prime numbers up to a given number N. The main idea of this algorithm is to sequentially remove all composite numbers that are divisible by 3, 5, 7, and so on. The numbers that remain unmarked are primes. But the Sieve of Eratosthenes is an algorithm, not an analytical formula. Although it is effective in practice for finding all prime numbers up to a given number N, it does not produce an exact analytical expression for the prime-counting function $\pi \left(N\right)$. The difficulty in obtaining an analytical formula is closely related to the complexity of implementing the sieving of numbers. Consequently, it becomes necessary to explore alternative computational or theoretical approaches in an effort to derive such a formula. In the following, the approach to approximate the calculation of $\pi \left(N\right)$ using analytic functions is presented.

Returning to the discussion of the Sieve of Eratosthenes algorithm, our objective is to quantify the number of composite numbers eliminated during the sieving process. To achieve this, we formulate an algorithm that computes the number of elements removed at each iteration, and after that we analyze the resulting pattern. Specifically, we examine how the number of remaining odds evolves throughout the procedure. Beyond a certain threshold, this amount stabilizes, indicating that only prime numbers remain. The number of these residual primes corresponds to the prime-counting function, denoted by the $\pi \left(N\right)$ function.

The pseudocode for this algorithm can be expressed as follows:

• 1) N = 5000
• 2) odds = list(3, N, step 2)
• 3) length_changes = [empty list]
• 4) i = 1
• 5) while (i < (some condition)):
•     (6) L = list(3∗(2∗i + 1), N, 2∗(2∗i) + 1)
•     (7) odds = difference(odds, L)
•     (8) length_changes.append(len(odds)(before)
•                     - len(odds)(after))
•     (9) i = i + 1

This is how the program can be constructed in a programming language.

• In line (1), a number is set — for example, N = 5000.
• In line (2), a list is created that contains all odd numbers from 3 to N.
• Then, on line (3), an empty list is initialized. This list will track how the length of the odd-number list changes as the composite numbers are removed.
• In block (5), a while loop runs at each step i (starting from 1), generating all the composite numbers divisible by $(2·i+1)$ from $3·(2·i+1)$ to N (line (6)), while a certain condition holds (the exact condition will be discussed in detail later).
• Lines (7) and (8) are crucial:

  -

  Line (7) performs a set difference between the list of remaining odd numbers and the composite numbers generated at the current step i. This effectively removes from the list all composite numbers divisible by $(2·i+1)$

  -

  Line (8) calculates how the length of the odd-numbers sequence changes, and this difference is stored in the list.

• Finally, line (9) proceeds to the next iteration i.

The array $odds$ after running this pseudocode contains all primes up to N, and its length is the total number of primes less than or equal to N, which is denoted as $\pi \left(N\right)$.

However, it is not necessary to find all the composite numbers for every step i up to N, because after a certain point no new composite numbers are generated, and the length of the odd-numbers array stops changing. Let this threshold be denoted as m. Then the condition in the while loop can be rewritten as $i=m$.

To determine m, we proceed with the following reasoning. The rows and columns in Figure 1 contain all numbers divisible by each odd number (each odd number can be generated with formula $2·i+1$ by letting i be equal to 1, 2, 3, etc.) within the range from 3 to N (we set $N=23$ in Figure 1 for simplicity).

• The first row contains all numbers divisible by 3, the second one by 5, etc.
• The table is symmetric with respect to the main diagonal, where the entries correspond to the squares of ($2·i+1$).
• All composite numbers - whose total count is to be determined - are located in one half of the matrix, along with diagonal elements that represent perfect squares.

Based on Figure 1, the number of composite numbers can be determined by simple reasoning:

• The number of values divisible by 3 is approximately ${\displaystyle \frac{N}{2·3}}$, considering only odd numbers.
• The count of values divisible by 5 is $\sim {\displaystyle \frac{N}{2·5}}-1$, since 15 has already been counted.
• The count of values divisible by 7 is $\sim {\displaystyle \frac{N}{2·7}}-2$, as 21 and 35 have already been included on the opposite side of the diagonal and should not be double-counted.

In general, for each successive odd number i, the number of previously counted values to subtract is $i-1$. The calculation continues until the subtraction term $i-1$ is equal to the number ${\displaystyle \frac{N}{2·(2·i+1)}}$, from which it is subtracted.

This simple reasoning can be expressed as follows:

$${\displaystyle \frac{N}{2·(2·i_m+1)}}=i_m-1$$

(13)

Where $i_m$ corresponds to the threshold m.

From condition (13) $i_m$ can be calculated as:

$$N=2·(2·i_m+1)(i_m-1)=2·(2·i_m^2-i_m-1)$$

(14)

$$(2·i_m^2-i_m-1)={\displaystyle \frac{N}{2}}$$

(15)

$$i_m^2\sim {\displaystyle \frac{N}{4}}$$

(16)

And finally:

$$i_m\sim {\displaystyle \frac{\sqrt{N}}{2}}$$

(17)

$$m\sim {\displaystyle \frac{\sqrt{N}}{2}}$$

(18)

We define m approximately as ${\displaystyle \frac{\sqrt{N}}{2}}$. The value m serves as a threshold beyond which the removal of composite numbers no longer contributes to the enumeration of composite numbers.

The pseudocode can then be rewritten as:

• 1) N = 5000
• 2) odds = list(3, N, step 2)
• 3) length_changes = [empty list]
• 4) i = 1
• ∗) m = N/(2(2i+1))
• 5) while (i < m):
•     (6) L = list(3∗(2∗i + 1), N, 2 ∗ (2∗i)+1)
•     (7) odds = difference(odds, L)
•     (8) length_changes.append(len(odds)(before)
•              - len(odds)(after))
•     (9) i = i + 1
•    (∗∗) m = N/(2(2i+1))

Two additional lines of code are introduced. The added line $(\ast )$ initializes the parameter m, which is subsequently reassigned to a new value (in another added line ($\ast \ast$)). Upon execution of the modified code, a sieving procedure can be applied to the list of odd numbers, effectively isolating the prime elements for an arbitrary N.

Figure 2 exhibits a notable oscillatory pattern and illustrates how the length of the odd number sequence evolves as the composite numbers are progressively removed. This behavior closely mirrors the characteristic form of the sinc function, defined as

$$\sin \mathrm{c}\left(x\right)={\displaystyle \frac{sin\left(x\right)}{x}}$$

(19)

To explore an approximation of the prime-counting function using the sinc function, we begin by considering the following definition:

$$\sin \mathrm{c}\left(wx\right)={\displaystyle \frac{sin\left(wx\right)}{wx}}$$

(20)

The function (20) exhibits oscillatory behavior and its frequency is governed by the parameter w. In the context of the sieving of composite numbers, specifically those divisible by 3 — we can naturally associate this frequency w with the reciprocal of the divisor, that is, $w={\displaystyle \frac{1}{3}}$. This natural choice reflects the periodicity of the removal of multiples of 3 in the sieve process. To appropriately scale the amplitude of the sinc function for our application, we introduce a normalization factor. Empirically, this factor can be determined approximately equal to ${\displaystyle \frac{1}{2}}·{\displaystyle \frac{N}{6}}$ (its determination is based on a direct comparison of variations of the length of the odds array depicted in a figure (2) with the behavior of the sinc function):

$$\mathrm{normalization}\mathrm{factor}={\displaystyle \frac{1}{2}}·{\displaystyle \frac{N}{6}}$$

(21)

Since the sinc function inherently oscillates between positive and negative values, and the sieving process involves the removal (not addition) of elements, we apply a modulus operation to ensure the non-negativity of the number of items removed. This transformation ensures that the function contributes only to the exclusion of composite numbers. Combining these components, the function can be rewritten as:

$$y\left(x\right)=|{\displaystyle \frac{1}{2}}·{\displaystyle \frac{N}{6}}·\sin \mathrm{c}\left({\displaystyle \frac{x}{3}}\right)|$$

(22)

The function $y\left(x\right)$ (22) should start from $x=1$. By definition, $\sin \mathrm{c}\left(0\right):=1$. However, in our case, $y\left(0\right)$ represents the number of the composite numbers divisible by 3, and this value should be subtracted initially. Therefore, formally we can define $y\left(0\right)={\displaystyle \frac{N}{3·2}}={\displaystyle \frac{N}{6}}$.

Plot a graph of the original function that calculates the variation in the length of an odd number sequence according to the pseudocode, along with the function $y\left(x\right)$ (for N = 5000) (see Figure 3 and the magnified image (Figure 4). Now we can compare the behavior of the original function (as we saw previously in Figure 2) and the function $y\left(x\right)$ (22). These two functions demonstrate very similar oscillation behavior; this leads one to believe that these functions are interconnected in some way. Note that upon magnifying the image, it becomes evident that the divergence between these functions grows with increasing i; however, we will defer a detailed study of this phenomenon to future work.

In case of $N=1000$, the function $y\left(x\right)$ fairly accurately follows the original function. If we compare the actual variations of odd-array’s length while removing the composite numbers and variations which are calculated by $y\left(x\right)$:

1) Original variations: [66, 37, 0, 20, 16, 0, 10, 8, 0, 6, 0, 0, 2, 1]

2) Calculated variations: [68, 34, 0, 17, 13, 0, 9, 8, 0, 6, 6, 0, 5, 4, 0]

Figure 5. The original function (blue line) of variations of odd numbers during of elimination of composite numbers (these variations are calculated according to the pseudocode and are hold in the list (lines 3 and 8 in the pseudocode)), and the function $y\left(x\right)$ (orange line) defined as (22) for $N=1000$).

Figure 5. The original function (blue line) of variations of odd numbers during of elimination of composite numbers (these variations are calculated according to the pseudocode and are hold in the list (lines 3 and 8 in the pseudocode)), and the function $y\left(x\right)$ (orange line) defined as (22) for $N=1000$).

The variations in the actual length of the odd-numbers array for $N=1000$ and the length calculated using the $y\left(x\right)$ function are quite similar. If we calculate the number of primes up to 1000 by eliminating composite numbers according to the original algorithm (i.e., by applying the direct sieving algorithm), and then calculate the number of primes using the $y\left(x\right)$ function for $N=1000$, we get the following results:

Calculation of primes by application of a direct sieving algorithm (as pseudocode describes):

$$p_1(N=1000)={\displaystyle \frac{N}{2}}-{\displaystyle \frac{N}{6}}-66-37-20-16-10-8-6-2-1=\mathbf{167}$$

Calculation of primes with function $y\left(x\right)$:

$$p_2(N=1000)={\displaystyle \frac{N}{2}}-{\displaystyle \frac{N}{6}}-68-34-17-13-9-8-6-5-4=\mathbf{169}$$

The actual number of primes is 167; the calculated number of primes with $y\left(x\right)$ function is 169.

The variations of odds-array’s length calculated by direct calculation and with $y\left(x\right)$ function are very similar. These findings imply that these two functions can exhibit correlated behavior or share underlying structural properties.

Now, rewrite the calculation the number of primes up to a given number N$\pi \left(N\right)$ using $y\left(x\right)$ function’s representation by a Formula (22)):

$$\pi \left(N\right)=N-{\displaystyle \frac{N}{2}}-{\displaystyle \frac{N}{6}}-\sum _{i=1}^m|{\displaystyle \frac{1}{2}}·{\displaystyle \frac{N}{6}}·\sin \mathrm{c}\left({\displaystyle \frac{i}{3}}\right)|$$

(23)

The subsequent analytical step involves considering that the original function, designed to evaluate changes in the length of the odd number array, yields zero at step i when $(2i+1)$ is not a prime number. Consequently, the function $y\left(x\right)$ should be redefined to include only those indices i for which $2·i+1=p$, where p is a prime. The summation process must therefore be restricted to prime-valued terms. However, given that elements divisible by 3 have already been excluded in a prior step, the primality condition must be verified for the subsequent value $2·(i+1)+1$, ensuring accurate filtration of non-prime contributions. This constraint must be incorporated into the formulation. Consequently, the revised expression should enforce the primality of ($2·i+1$) as a necessary condition for inclusion in the summation.

$$\pi \left(N\right)=N-{\displaystyle \frac{N}{2}}-{\displaystyle \frac{N}{6}}-\sum _{i_p}^m|{\displaystyle \frac{1}{2}}·{\displaystyle \frac{N}{6}}·\sin \mathrm{c}\left({\displaystyle \frac{i_p}{3}}\right)|=$$

$$=N-{\displaystyle \frac{N}{2}}-{\displaystyle \frac{N}{6}}-{\displaystyle \frac{1}{2}}·{\displaystyle \frac{N}{6}}·\sum _{i_p}^m\left|\sin \mathrm{c}\left({\displaystyle \frac{i_p}{3}}\right)\right|$$

(24)

where $i_p={\displaystyle \frac{p-1}{2}}$, indicating that the summation is restricted to prime numbers only; the upper bound of the summation is defined by $m={\displaystyle \frac{\sqrt{N}}{2}}$, which specifies the limiting value to which the summation is carried out.

The Formula (24) is an important one. First, we subtract all even numbers from N (their total value is ${\displaystyle \frac{N}{2}}$), then all numbers that are divisible by 3 (their total value is ∼${\displaystyle \frac{N}{2·3}}={\displaystyle \frac{N}{6}}$), and finally we perform a sequential subtraction of all values $y\left(i\right)$, with the condition that the number $2·(i+1)+1$ should be prime.

However, the function (24) yields an approximate solution of $\pi \left(N\right)$ rather than an exact value. To reflect this, $\pi \left(N\right)$ in the Formula (24) can be reformulated as follows:

$$\pi \left(N\right)=F\left(N\right)+Error\left(N\right)$$

(25)

In Formula (25), $F\left(N\right)$ can be interpreted as a simple but dominant function that is very close to the real value of $\pi \left(N\right)$. Based on this, we can decompose the function $\pi \left(N\right)$ defined in (24) into two constituent components:

$$\pi \left(N\right)\sim N-{\displaystyle \frac{N}{2}}-{\displaystyle \frac{N}{6}}-{\displaystyle \frac{1}{2}}·{\displaystyle \frac{N}{6}}·\sum _{i_p}^m\left|\sin \mathrm{c}\left({\displaystyle \frac{i_p}{3}}\right)\right|\sim$$

$$\sim N-{\displaystyle \frac{N}{2}}-{\displaystyle \frac{N}{6}}-{\displaystyle \frac{1}{2}}·{\displaystyle \frac{N}{6}}·\sum _{i_p}^m\left|\right(f\left(i_p\right)+error\left(i_p\right)\left)\right|$$

(26)

Let $f\left(i_p\right)$ denote the dominant component that closely approximates the variations in the length of the array of odd numbers, and let $error\left(i_p\right)$ represent an oscillatory error value corresponding to the variation at the step ($i_p$). At this stage, we are not concerned with the explicit analytical forms of $f\left(i_p\right)$ and $error\left(i_p\right)$; rather, we assume the existence of such functions that capture the described behavior.

Decompose the function $\pi \left(N\right)$ into two distinct components:

$$\pi \left(N\right)=N-{\displaystyle \frac{N}{2}}-{\displaystyle \frac{N}{6}}-{\displaystyle \frac{1}{2}}·{\displaystyle \frac{N}{6}}·\sum _{i_p}^m\left|\right(f\left(i_p\right)+error\left(i_p\right)\left)\right|=$$

$$=N-{\displaystyle \frac{N}{2}}-{\displaystyle \frac{N}{6}}-{\displaystyle \frac{1}{2}}·{\displaystyle \frac{N}{6}}·\sum _{i_p}^{m}f\left(i_p\right)-{\displaystyle \frac{1}{2}}·{\displaystyle \frac{N}{6}}·\sum _{i_p}^{m}error\left(i_p\right)$$

(27)

The term $N-{\displaystyle \frac{N}{2}}-{\displaystyle \frac{N}{6}}-{\displaystyle \frac{1}{2}}·{\displaystyle \frac{N}{6}}·{\sum }_{i_p}^{m}f\left(i_p\right)$ represents the dominant component which we denote as $F\left(N\right)$, and the term ${\displaystyle \frac{1}{2}}·{\displaystyle \frac{N}{6}}·{\sum }_{i_p}^{m}error\left(i_p\right)$ represents the error component responsible for capturing the error behavior, which we denote as $E\left(N\right)$:

$$F\left(N\right)=N-{\displaystyle \frac{N}{2}}-{\displaystyle \frac{N}{6}}-{\displaystyle \frac{1}{2}}·{\displaystyle \frac{N}{6}}·\sum _{i_p}^{m}f\left(i_p\right)$$

(28)

$$E\left(N\right)={\displaystyle \frac{1}{2}}·{\displaystyle \frac{N}{6}}·\sum _{i_p}^{m}error\left(i_p\right)$$

(29)

Since summation in (29) is restricted to prime numbers, we utilize the Prime Number Theorem (PNT) to estimate a number of terms in the ${\sum }_{i_p}^{m}error\left(i_p\right)$ sequence in $E\left(N\right)$. Specifically, the number of primes up to m can be approximated by $\pi \left(m\right)={\displaystyle \frac{m}{lnm}}$. This approximation allows us to estimate the number of terms involved in the summation over primes up to the limit m:

$$\pi \left(m\right)=\pi \left({\displaystyle \frac{\sqrt{N}}{2}}\right)={\displaystyle \frac{\frac{\sqrt{N}}{2}}{ln\frac{\sqrt{N}}{2}}}={\displaystyle \frac{\sqrt{N}}{2(ln\sqrt{N}-ln2)}}\sim {\displaystyle \frac{\sqrt{N}}{lnN}}$$

(30)

And now rewrite $\pi \left(N\right)$ as:

$$\pi \left(N\right)=F\left(N\right)-E\left(N\right)=F\left(N\right)-{\displaystyle \frac{1}{2}}·{\displaystyle \frac{N}{6}}·\sum _{i_p}^{m}error\left(i_p\right)\sim$$

$$\sim F\left(N\right)-{\displaystyle \frac{\sqrt{N}}{lnN}}·E_{mean}\left(N\right)$$

(31)

where $E_{mean}\left(N\right)$ is a function of N that accounts for a mean error across all $error\left(i_p\right)$ and ${\displaystyle \frac{\sqrt{N}}{lnN}}$ is an estimation of a total number of terms in $E\left(N\right)$.

Finally, rewrite $\pi \left(N\right)$ as

Formula (32) implies that the prime-counting function $\pi \left(N\right)$ can be expressed in two terms: a function $F\left(N\right)$, which closely approximates $\pi \left(N\right)$ - and an average error component ${\displaystyle \frac{\sqrt{N}}{lnN}}·E_{mean}\left(N\right)$, representing the cumulative deviation from the expected distribution. Furthermore, the error term is expected to exhibit an oscillatory behavior, analogous to that of the sinc function, reflecting the inherent fluctuations in the distribution of prime numbers. For the purposes of this analysis, we shall abstract away the intrinsic nature of the function $F\left(N\right)$. It is essential that $F\left(N\right)$ represents the dominant contribution to $\pi \left(N\right)$ and $E_{mean}\left(N\right)$ represents the error term.

Formula (32) is the most significant implication derived from the preceding analysis. In fact, it allows us to conclude that the error term in the prime-counting function $\pi \left(N\right)$ should have a behavior similar to ${\displaystyle \frac{\sqrt{N}}{lnN}}·E_{mean}\left(N\right))$ regardless of the $F\left(N\right)$ function, where the $E_{mean}\left(N\right)$ function represents an oscillatory error term with bound coefficient ${\displaystyle \frac{\sqrt{N}}{lnN}}$.

Return to the Chebyshev prime-counting function $\psi \left(x\right)$, which was described in the introduction (11). 35 years after Riemann published his paper, the German mathematician Hans von Mangoldt studied $\psi \left(x\right)$ and found a relatively simple explicit formula for it [5]:

$$\psi \left(x\right)=x-log\left(2\pi \right)-{\displaystyle \frac{1}{2}}log(1-{\displaystyle \frac{1}{x^2}})-\sum _{\rho }{\displaystyle \frac{x^{\rho }}{\rho }}$$

(33)

This formula consists of four terms: the first one is simply x, the second one $log2\pi$ is a constant, the third one ${\displaystyle \frac{1}{2}}log(1-{\displaystyle \frac{1}{x^2}})$, which involves the sum of trivial zeros of the zeta function and becomes negligible for large value of x, and the fourth and the important one is the sum of ${\displaystyle \frac{x^{\rho }}{\rho }}$ over all the non-trivial zeros of the zeta function.

Divide equation (12) by $logx$:

$$|{\displaystyle \frac{\psi \left(x\right)}{logx}}-{\displaystyle \frac{x}{logx}}|<{\displaystyle \frac{\sqrt{x}logx}{8\pi }}$$

(34)

Comparison of (34) and (11) gives the following:

$$|{\displaystyle \frac{\psi \left(x\right)}{logx}}-{\displaystyle \frac{x}{logx}}|\sim |\pi \left(x\right)-\mathrm{li}\left(\mathrm{x}\right)|$$

(35)

From (33) after exclusion of $log2\pi$ and ${\displaystyle \frac{1}{2}}log(1-{\displaystyle \frac{1}{x^2}})$ an expression for $\psi \left(x\right)$ can be rewritten as a relatively simple explicit formula:

$$\psi \left(x\right)\sim x-\sum _{\rho }{\displaystyle \frac{x^{\rho }}{\rho }}$$

(36)

Divide (36) by $logx$:

$${\displaystyle \frac{\psi \left(x\right)}{logx}}\sim {\displaystyle \frac{x}{logx}}-{\displaystyle \frac{1}{logx}}\sum _{\rho }{\displaystyle \frac{x^{\rho }}{\rho }}$$

(37)

$${\displaystyle \frac{\psi \left(x\right)}{logx}}-{\displaystyle \frac{x}{logx}}\sim -{\displaystyle \frac{1}{logx}}\sum _{\rho }{\displaystyle \frac{x^{\rho }}{\rho }}$$

(38)

From (35) derive an expression for the difference ($\pi \left(x\right)-\mathrm{li}\left(\mathrm{x}\right)$):

$${\displaystyle \frac{\psi \left(x\right)}{logx}}-{\displaystyle \frac{x}{logx}}=\pi \left(x\right)-\mathrm{li}\left(\mathrm{x}\right)\sim -{\displaystyle \frac{1}{logx}}\sum _{\rho }{\displaystyle \frac{x^{\rho }}{\rho }}$$

(39)

$$\pi \left(x\right)\sim \mathrm{li}\left(\mathrm{x}\right)-{\displaystyle \frac{1}{logx}}\sum _{\rho }{\displaystyle \frac{x^{\rho }}{\rho }}$$

(40)

The expression (40) shows that $\pi \left(x\right)$ is estimated as $\mathrm{li}\left(x\right)$ minus the sum of all non-trivial zeros of the zeta function divided by $logx$. If RH is true, all non-trivial zeros have a real part that is ${\displaystyle \frac{1}{2}}$, and thus each non-trivial zero $\rho$ can be written as $\rho ={\displaystyle \frac{1}{2}}+i\sigma$.

Using this, rewrite the sum ${\sum }_{\rho }{\displaystyle \frac{x^{\rho }}{\rho }}$:

$$\sum _{\rho }{\displaystyle \frac{x^{\rho }}{\rho }}=\sum _{\rho }{\displaystyle \frac{x^{\frac{1}{2}+i\sigma }}{\frac{1}{2}+i\sigma }}=\sum _{\rho }{\displaystyle \frac{\sqrt{x}·x^{i\sigma }}{\frac{1}{2}+i\sigma }}=\sqrt{x}\sum _{\rho }{\displaystyle \frac{x^{i\sigma }}{\frac{1}{2}+i\sigma }}=\sqrt{x}·g\left(x\right)$$

(41)

The $g\left(x\right)$ function in (41) is a sum of imaginary parts of the non-trivial zeros of the zeta function.

Based on (41), the expression for the sum ${\sum }_{\rho }{\displaystyle \frac{x^{\rho }}{\rho }}$ can be written in the simplified version:

$$\sum _{\rho }{\displaystyle \frac{x^{\rho }}{\rho }}=\sqrt{x}·g\left(x\right)$$

(42)

And finally, Formula (40) can be written as:

$$\pi \left(x\right)\sim \mathrm{li}\left(\mathrm{x}\right)-{\displaystyle \frac{1}{logx}}\sum _{\rho }{\displaystyle \frac{x^{\rho }}{\rho }}=\mathrm{li}\left(\mathrm{x}\right)-{\displaystyle \frac{\sqrt{x}}{logx}}g\left(x\right)$$

(43)

where $g\left(x\right)$ is equal to the sum of imaginary parts of non-trivial zeros of the zeta function ${\sum }_{\rho }{\displaystyle \frac{x^{\rho }}{\rho }}$.

Now we have two formulas for $\pi \left(x\right)$:

The formula (45) is derived from the analysis of oscillations of the length of odd-numbers array, where $E_{mean}\left(N\right)$ is an oscillation function and $F\left(x\right)$ is a function that is close to the true value of $\pi \left(x\right)$ and can be considered as an analog for li(x) (because function li(x) is the closest known approximation of the prime-counting function). More interesting is a $g\left(x\right)$ term which aggregates a sum of imaginary parts of non-trivial zeros.

If we rewrite $x^{i\sigma }$ as:

$$x^{i\sigma }=e^{ln{x}^{i\sigma }}=e^{i\sigma lnx}$$

(46)

and then the sum of non-trivial zeros can be rewritten as next:

$$\sum _{\rho }{\displaystyle \frac{x^{\rho }}{\rho }}=\sum _{\rho }{\displaystyle \frac{x^{\frac{1}{2}}{x}^{i\sigma }}{\frac{1}{2}+i\sigma }}\sim \sum _{\sigma }{\displaystyle \frac{\sqrt{x}{e}^{i\sigma lnx}}{i\sigma }}=\sqrt{x}\sum _{\sigma }{\displaystyle \frac{e^{i\sigma lnx}}{i\sigma }}$$

(47)

The function $e^{i\sigma lnx}$ also has an oscillating behavior, because it can be represented by Euler’s formula:

$$e^{i\sigma lnx}=cos(i\sigma lnx)+isin(i\sigma lnx)$$

(48)

Rewrite the sum (47) by using (48):

$$\sum _{\sigma }{\displaystyle \frac{e^{i\sigma lnx}}{i\sigma }}=\sum _{\sigma }{\displaystyle \frac{cos(i\sigma lnx)+isin(i\sigma lnx)}{i\sigma }}=$$

$$=\sum _{\sigma }{\displaystyle \frac{cos(i\sigma lnx)}{i\sigma }}+\sum _{\sigma }{\displaystyle \frac{sin(i\sigma lnx)}{\sigma }}$$

(49)

Riemann’s hypothesis expresses the error in the prime number theorem as the sum of waves with wavelength determined by the locations of the zeros of $\zeta \left(s\right)$. In [7] the following approximate expression is given to demonstrate this:

$${\displaystyle \frac{\mathrm{Li}\left(\mathrm{x}\right)-\pi \left(x\right)}{\sqrt{x}/lnx}}\sim 1+2\sum _{\sigma }{\displaystyle \frac{sin(\sigma lnx)}{\sigma }}$$

(50)

The zeros of the Riemann zeta function are placed so precisely that they anticipate the locations of the primes. Moreover, they do it through the Fourier series type term in (50) in which the values of the zeros determine both the amplitude and the wavelength - analog of the sine-log waves. If the Riemann hypothesis is false, the Formula (50) would be marred with more complicated coefficients for the sines that are functions of x, which are different from $\sqrt{x}/lnx$.

Regrouping the terms, we can rewrite our formula (45), which was previously derived from the analysis of variations of the odd number quantity during the removal of composite numbers, in the following form of (50):

$${\displaystyle \frac{F\left(x\right)-\pi \left(x\right)}{\sqrt{x}/lnx}}\sim E_{mean}\left(N\right)$$

(51)

The function $E_{mean}\left(N\right)$ is the mean of errors up to the limit m (see Formula (31)) and also has sinusoidal oscillation behavior in the form of functions sinc, as we saw previously (26). Based on the similarity of (50) and (51), we can propose that this $E_{mean}\left(N\right)$ function describes a sum of imaginary parts of non-trivial zeros of the zeta function. The more interesting for us is that this formulation concludes that the coefficient in our Formula (51) has a value of ${\displaystyle \frac{\sqrt{x}}{logx}}$, which we have derived from an estimation of the limit to count the composite numbers, as we established earlier. But, in formulas (47) and (50), this coefficient is derived from the proposal that all non-trivial zeros have a real part equal to ${\displaystyle \frac{1}{2}}$, as proposed by the Riemann Hypothesis. Moreover, the similarity of formulas (51), (50) and (44) may also provide supporting evidence for the validity of the RH.

3. Conclusions

In this work, an algorithm for calculating the prime counting function $\pi \left(x\right)$ is proposed by determining the number of odd numbers remaining after removing the composite ones. This method generalizes the classical Sieve of Eratosthenes algorithm by not only identifying and excluding composite numbers but also quantifying the number of such eliminations during the sieving process. By applying this approach, it was found that these variations follow an oscillation pattern governed by the sinc function. The further analysis of such oscillations allows us to compose the function $\pi \left(x\right)$ by two terms: the main term, which makes the largest contribution to the distribution of the primes, and the error term, which is responsible for the accumulation of errors during the calculation of the functions sinc across the prime numbers p up to the limit m defined as $\sqrt{x}/2$. This leads to the conclusion that the bound coefficient for the error term in PNT should be equal to $\sqrt{x}/logx$, which correlates with an estimation of this coefficient derived under the assumption of the truth of the Riemann hypothesis. We hope that this perspective stimulates future work toward formalizing this approach to uncover deeper connections to the zeta function and prime number theory at large scale.