Saturday, September 12, 2015

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

No comments:

Post a Comment