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
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
No comments:
Post a Comment