Sunday, September 13, 2015

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

No comments:

Post a Comment