Counting primes - Prime Counting function - Part III

In the lines of Divisor Summatory function, let us define a Prime Summatory function $\pi_k(x)$. That is,

$\pi_k(x)=\displaystyle\sum_{p\leq x}p^k$

where $p$ is prime.

Then we can use the same function we used in Counting primes - Prime Counting function - Part II to calculate $\pi_k(x)$. We can define a function $\varphi_k(x,a)$ which gives the sum of the $k$th powers of integers between $2$ and $x$ (both inclusive), that remains unmarked after sieving with the first $a$ primes. Then as before,

$\pi_k(x) = \varphi_k(x,\pi(\sqrt{x}))$

and

$\varphi_k(x,a)=
\begin{cases}
\varphi_k(x,a-1)-p_a^k(\varphi_k(x/p_a,a-1)-\varphi_k(p_a-1,a-1)),&\text{if }x\geq p_a^2\\
\varphi_k(x,\pi(\sqrt{x})),&\text{otherwise}
\end{cases}$

such that $\varphi_k(0,a)=\varphi_k(1,a)=0$ and $\varphi_k(x,0)=2^k + 3^k + \cdots + x^k$.

Even the same Mathematica code used in the previous post with a little modification does the job. For example, the following Mathematica code gives the sum of primes less than a given value.

Clear[S, PrimeP, Su]
S[n_, m_] :=
 S[n, m] =
  Which[n <= 1, 0, m == 0, n - 1, n < Power[Prime[m], 2], PrimeP[n],
   True, S[n,
     m - 1] - (S[Floor[n/Prime[m]], m - 1] - S[Prime[m] - 1, m - 1])]
PrimeP[n_] :=
 PrimeP[n] = If[n <= 3, n - 1, S[n, PrimeP[Floor[Sqrt[n]]]]]
Su[n_, m_] :=
 Su[n, m] =
  Which[n <= 1, 0, m == 0, n (n + 1)/2 - 1, n < Power[Prime[m], 2],
   Su[n, PrimeP[Floor[Sqrt[n]]]], True,
   Su[n, m - 1] -
    Prime[m] (Su[Floor[n/Prime[m]], m - 1] - Su[Prime[m] - 1, m - 1])]
n = Power[10, 9];
Timing[Block[{RecursionLimit = Power[10, 8]}, Su[n, PrimeP[n]]]]

This code computes the value of $\pi_1(10^8)$ in $7$ seconds and the value of $\pi_1(10^9)$ in about $30$ seconds. Though the code does not have any significant changes from that of $\pi(x)$, the code $\pi_1(x)$ still relies on $\pi(x)$ which increases the calculation time.

Won't it be nice if we can somehow eliminate that and thereby increase the performance? We'll see about that in the next post.


Until then
Yours Aye
Me

Counting primes - Prime Counting function - Part II

We discussed the Legendre's formula associated with the prime counting function in the previous post. This post discusses an alternative formulation whose importance will be apparent when we discuss the SumOfPrimes function.

Define a new function $\varphi(x,a)$ such that it counts the number of integers between $2$ and $x$ that remain unmarked after sieving with the first $a$ primes. Such a function will be related to the prime-counting function as

$\pi(x) = \varphi(x,\pi(\sqrt{x}))$

This function satisfies the recurrence relation

$\varphi(x,a)=
\begin{cases}
\varphi(x,a-1)-(\varphi(x/p_a,a-1)-\varphi(p_a-1,a-1)),&\text{if }x\geq p_a^2\\
\pi(x),&\text{otherwise}
\end{cases}$

such that $\varphi(0,a)=\varphi(1,a)=0$ and $\varphi(x,0)=x-1$.

Note that whereas $\phi(10,2)=3$ (the three numbers being $1, 5$ and $7$), $\varphi(10,2)=4$ (the four numbers are $2, 3, 5$ and $7$). Because when sieve with the $2$, every even number except $2$ goes marked, with $3$ every mulitple of $3$ except $3$ itself goes marked and so on.

The following Mathematica code uses these properties of $\varphi(x,a)$ to calculate $\pi(x)$.

Clear[S, PrimeP, Su]
S[n_, m_] :=
 S[n, m] =
  Which[n <= 1, 0, m == 0, n - 1, n < Power[Prime[m], 2], PrimeP[n],
   True, S[n,
     m - 1] - (S[Floor[n/Prime[m]], m - 1] - S[Prime[m] - 1, m - 1])]
PrimeP[n_] :=
 PrimeP[n] = If[n <= 3, n - 1, S[n, PrimeP[Floor[Sqrt[n]]]]]
Su[n_, m_] :=
 Su[n, m] =
  Which[n <= 1, 0, m == 0, n (n + 1)/2 - 1, n < Power[Prime[m], 2],
   Su[n, PrimeP[Floor[Sqrt[n]]]], True,
   Su[n, m - 1] -
    Prime[m] (Su[Floor[n/Prime[m]], m - 1] - Su[Prime[m] - 1, m - 1])]
n = Power[10, 9];
Timing[Block[{RecursionLimit = Power[10, 8]}, PrimeP[n]]]

This code computes $\pi(10^7)$ in a sec, takes approx. 4 secs for $\pi(10^8)$ and 18 secs for $\pi(10^9)$. As we can see, there is no significant increase in performance between the two codes. But as I said at the beginning of this post, this function can be easily tweaked to find the SumOfPrimes function because of its lack of dependence on $a$ in the second case of the recurrence.

We'll see about that in detail in the next post.


Until then
Yours Aye
Me

Counting primes - The Prime-counting function - Legendre's Phi function

The prime-counting function, denoted as $\pi(x)$, counts the number of primes less than or equal to some value x. Though there are many formulas to calculate the value of this function approximately, finding the exact value is a challenging problem, especially for beginner level programmer.

Apart from naive counting, the next best formula used for counting primes was given by Legendre. You can read about it in the wiki page. This formula introduces a new function called the Legendre's phi function, which is then used to compute $\pi(x)$.

However, in this post we'll use a closely related function that is slightly more efficient than the Lengendre's formula given in Wiki or Math-world page.

According to Math-world, Legendre's Phi function, denoted as $\phi(x,a)$, counts the number of integers less or equal to $x$ which are not divisible by any of the first $a$ primes. This function is then related to the prime-counting function as

$\phi(x,\pi(\sqrt{x}))=\pi(x)-\pi(\sqrt{x})+1$

Also, it follows the following recurrence relation

$\phi(x,a)=
\begin{cases}
\phi(x, a-1) - \phi(x/p_a, a-1), &\text{if }x\geq p_a^2 \\
\pi(x)-a+1,&\text{otherwise}
\end{cases}$

where $p_a$ denotes the $a$th prime. Also, $\phi(0, a) = 0$, $\phi(1, a) = 1$ and $\phi(x, 0)=x$.

The following Mathematica code uses these relations to compute $\pi(x)$.

Clear[S, PrimeP]
S[n_, m_] :=
 S[n, m] =
  Which[n <= 1, n, m == 0, n, n < Power[Prime[m], 2],
   PrimeP[n] - m + 1, True, S[n, m - 1] - S[Floor[n/Prime[m]], m - 1]]
PrimeP[n_] :=
 PrimeP[n] =
  If[n <= 3, n - 1,
   S[n, PrimeP[Floor[Sqrt[n]]]] + PrimeP[Floor[Sqrt[n]]] - 1]
n = Power[10, 9];
Timing[Block[{RecursionLimit = Power[10, 8]}, PrimeP[n]]]

The code takes approx a second to calculate $\pi(10^7)$, 4 secs for $\pi(10^8)$ and 20 secs for $\pi(10^9)$. Though these timings are reasonable for Mathematica, and will speed up considerably if we do in C, we'll see an alternative formulation in the next post.

Until then
Yours Aye
Me