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

C-Code for Merten's Function

The C-Code for Merten's function proceeds much in the same way as that of Totient Summatory function and uses the formula derived in Merten's function. This code evaluates the value of $M(10^{11})$ in approx. $8$ seconds and $M(10^{12})$ in approx. $80$ seconds. Here is the code.

#include<stdio.h>
#include<math.h>

#define M 60000000
#define Mod 1000000007

int MertCache[M+1], Mert[M+1];
char Mu[M+1], flag[M+1];

void precalc()
{
    int i, j;
    Mu[1] = 1;
    Mert[1] = 1;
    MertCache[1] = 56789;
    for(i = 2; i <= M; i++) {
        if(!flag[i])
            for(j = i; j <= M; j += i) {
                if( !flag[j] )
                    Mu[j] = 1;
                if(j/i%i == 0)
                    Mu[j] = 0;
                else
                    Mu[j] = -Mu[j];
                flag[j] = 1;
        }
        Mert[i] = Mert[i-1] + Mu[i];
        MertCache[i] = 56789;
    }
}

long long Merten(long long n, long long k)
{
    long long m = n/k;
    if(m <= M)
        return Mert[m];
    else if(MertCache[k] != 56789)
        return MertCache[k];
    long long p = sqrt(m), q = m/(p+1), i, res = 0;
    for(i = 2; i <= p; i++)
        res += Merten(n, k*i);
    for(i = 1; i <= q; i++)
        res = res + (m/i) * Mu[i];
    res = 1 - res;
    res += p * Merten(n, k*(p+1));
    return MertCache[k] = res;
}

int main()
{
    precalc();
    printf("%lld", Merten(100000000000LL,1));
    return 0;
}

Note: $\text{Merten}(n,k) = M\left(\left\lfloor\dfrac{n}{k}\right\rfloor\right)$ in the code.

This code doesn't include the Phi function at all. But if you plan to use both the Totient summatory function and the Merten's function in the same code, I think code may slow down. I'vnt tried it yet but I'll also post the code that uses both the functions.

UPDATE (13 May 2017): As in Totient Summatory function, I made an iterative version of the Merten's function. You can check the performance of the code in this Online IDE.

from timeit import default_timer

start = default_timer()

n = 10 ** 6
s = int(n ** 0.5)


def sieve():
    Lows, Flag = [0] + [1] * s, [0] * (1 + s)
    Lows[1] = 1
    for i in range(2, 1 + s):
        if Flag[i] == 0:
            for j in range(i, 1 + s, i):
                Flag[j] = 1
                if (j // i % i) == 0:
                    Lows[j] = 0
                else:
                    Lows[j] *= -1
        Lows[i] += Lows[i - 1]
    return Lows


Lows = sieve()
Highs = [0] * (1 + s)

for k in range(s, 0, -1):
    tempn = n // k
    u = int(tempn ** 0.5)
    v = tempn // (1 + u)
    res = 1
    for i in range(1, 1 + v):
        res -= (Lows[i] - Lows[i - 1]) * (tempn // i)
    for i in range(2, 1 + u):
        res -= Highs[k * i] if k * i <= s else Lows[tempn // i]
    i = 1 + u
    res += Lows[tempn // i] * u
    Highs[k] = res

print(Highs[1])
print(default_timer() - start)

Note: $\text{Highs}[k] = M\left(\left\lfloor\dfrac{n}{k}\right\rfloor\right)$ in the code. Also, the value $-962982$ is purely arbitrary just to make sure that value has not been touched yet.


Until then
Yours Aye
Me

C Code for Totient Summatory function

I thought of creating a post on Dirichlet Generating functions which would be right continuation of the direction my blog is heading. But this post is about the C code for Totient Summatory function.

The code uses the algorithm described in Totient Summatory function. If you haven't visited that post, I suggest you to do so. The code computes the answer for $10^{11}$ in approx. $7$ seconds and for $10^{12}$ in approx. $25$ seconds. Not bad compared to Mathematica. Here's the code for you.

#include<stdio.h>
#include<math.h>

#define M 60000000
#define Mod 1000000007

int Phi[M+1], SumPhi[M+1], TotCache[M+1];

void precalc()
{
    int i, j;
    Phi[1] = 1;
    SumPhi[1] = 1;
    for(i = 2; i <= M; i++) {
        if(!Phi[i])
            for(j = i; j <= M; j += i) {
                if(!Phi[j])
                    Phi[j] = j;
                Phi[j] = (Phi[j] / i) * (i - 1); }
        SumPhi[i] = (SumPhi[i-1] + Phi[i]) % Mod; }
}

long long SumNat(long long n)
{
    long long p = n, q = n+1;
    if(p%2)
        q /= 2;
    else
        p /= 2;
    return (p % Mod) * (q % Mod) % Mod;
}

long long TotientSum(long long n, long long k)
{
    long long m = n / k;
    if(m <= M)
        return SumPhi[m];
    else if(TotCache[k] > 0)
        return TotCache[k];
    long long p = sqrt(m), q = m / (p+1), i, res = 0;
    for(i = 2; i <= p; i++) {
        res += TotientSum(n, k * i);
        if(res >= Mod)
            res %= Mod;
    }
    for(i = 1; i <= q; i++) {
        res = res + ((m / i) % Mod * Phi[i]) % Mod;
        if(res >= Mod)
            res %= Mod;
    }
    res = SumNat(m) - res;
    res += (p % Mod) * TotientSum(n, k * (p + 1));
    res %= Mod;
    TotCache[k] = res;
    return res;
}

int main()
{
    precalc();
    printf("%lld",TotientSum(1000000000000LL,1));
    return 0;
}

Note: $\text{TotientSum}(n,k)=\Phi\left(\left\lfloor\dfrac{n}{k}\right\rfloor\right)$ in the code.

Just posting a code wouldn't do much good either for the blogger or the reader I think from the experience I've had myself. But am not gonna no through the complete code. Just a few details which we should be paying attention to.

First, the 'precalc' function. It's just a word-by-word translation of the algorithm that we saw on Optimized Precomputation.

Second is the 'SumNat' function. The important thing about this function is that we ought to keep the values from overflowing. For example, I first did $\text{(p%Mod * q%Mod)%Mod}$ which was perfectly okay looking except that it overflows.

Third, the 'TotientSum' function, is pure memoization. It'll take a separate post to completely explain it. Maybe I'll do that later. But for now I think this should suffice. I'll try to come up with the C-code for Merten's function in the next post.

Last is the limit you choose to sieve. That's the real trade-off you have to make between space and time. If you choose to sieve more, it's gonna take more space but the final result will be quick, whereas if you choose to memoize more you consume less space at cost of being time expensive. I made some trial and error to arrive at the above sieve limit which kinda minimized the computation time.

UPDATE(13 Mar 2017): I made an iterative version of the Totient summartory function in Python. You can check it out in this Online IDE.

from timeit import default_timer

start = default_timer()

n = 1000000000
s = int(n ** 0.5)
LowPhis = [k for k in range(1 + s)]
HiPhis = [0] * (1 + s)

for p in range(2, 1 + s):
    if p == LowPhis[p]:
        for k in range(p, 1 + s, p):
            LowPhis[k] = LowPhis[k] - LowPhis[k] // p
    LowPhis[p] += LowPhis[p - 1]

for k in range(s, 0, -1):
    tempn = n // k
    res = tempn * (1 + tempn) // 2
    u = int(tempn ** 0.5)
    v = tempn // (1 + u)
    for i in range(1, 1 + v):
        res -= (tempn // i) * (LowPhis[i] - LowPhis[i - 1])
    for i in range(2, 1 + u):
        res -= (HiPhis[i * k] if i * k <= s else LowPhis[tempn // i])
    i = 1 + u
    res += u * LowPhis[tempn // i]
    HiPhis[k] = res

print(HiPhis[1])
print(default_timer() - start)

Note: $\text{HiPhis}[k]=\Phi\left(\left\lfloor\dfrac{n}{k}\right\rfloor\right)$ in the code.

I like iterative versions in that they don't suffer from the annoying maximum-recursion-depth error. I think this is a little slower than the recursive version but looks well nonetheless in Python.


Until then
Yours Aye
Me

An Optimized Pre Computation

This post is an extension to make an optimization in precomputing the values of Totient and Mobius function. So make sure you first go through the previous post.

The extra initialization required for Algorithm 2 made me think for a while to eliminate the 'initializing' for loop. So I made (or atleast tried to make) a little optimization there. I realized that we don't need that loop. Instead, I initialized the values in the for loop only when the EulerPhi was zero. In other words, when we first reach a number, we do the initialization which should be sufficient enough. The following algorithm do that.

Algorithm 3

function PRECOMPUTATION(n)

            flag  $\gets$ [0] * (n)                                                                                            

            EulerPhi $\gets$ [0] * (n)

            Mu $\gets$ [0] * (n)

            Primes $\gets$ [0] * (n)

            pos $\gets$ 1

            for i = 2 to n do

if flag[i] = 0 then

                                    for j = i to n do

                                                flag[j] = 1

                                                if EulerPhi[j] = 0 then

                                                            EulerPhi[j] $\gets$ j

                                                            Mu[j] $\gets$ 1

                                                endif

                                                EulerPhi[j] $\gets$ ( EulerPhi[j] / i ) * ( i– 1 )

                                                if j/i%i = 0 then

                                                            Mu[j] $\gets$ 0

                                                else

                                                            Mu[j] $\gets$ Mu[j] * (-1)

                                                end if

                                                j $\gets$ j + i

                                    end for

                                    Primes[pos] $\gets$ i

                                    pos $\gets$ pos +1

                           end if

           end for

          return Primes

          return EulerPhi

          return Mu

end function

This pretty much completes our discussion of Precomputations. These algorithms can be effectively converted into codes in a programming language of your choice.

Yours Aye

Me

Pre Computation

The pre-computation of Totient function and Mobius function is one of the prerequisites of all the algorithms that we have seen so far. In fact, this pre-computation would form an excellent exercise in understand the standard Sieve of Erasthones.  In our sieving, we'll not only find the primes, but also the two aforementioned functions.

Let's first see the algorithm for Erasthones' sieve. It's not difficult understand or code. From there we'll build our Totient and Mobius functions.

Algorithm 1

function PRIMES(n)

            flag  $\gets$ [0] * (n)

            Primes $\gets$ [0] * (n)

            pos $\gets$ 1

            for i = 2 to n do

if flag[i] = 0 then

                                    for j = i to n do

                                                flag[j] $\gets$ 1
                                                j $\gets$ j + i
                                    end for

                                    Primes[pos] $\gets$ i

                                    pos $\gets$ pos +1

                           end if

            end for

            return Primes

end function

As you can see here, we start our algorithm by initializing a flag of 0's. Starting from 2, whenever we first encounter a 0, we label it as a prime. We then mark it's multiples as composites by changing their flags to 1. We then move on to find the next 0 and repeat the process.

Well, in order to bring in the totient function, we first note that the totient function of a particular number $n$ is set to its own value and changes only by its prime factorization. We encounter all the multiples of a prime in the inner for loop and use this loop to change the totient value as appropriate.

The logic applies to Mobius function as well. We initially set the value to be 1 and multiply it by -1 whenever we reach a number by a different prime. In addition, for Mobius we also have to check if its divisible by a square, and if it is, we change it to 0.

Algorithm 2

function PRECOMPUTATION(n)

            flag  $\gets$ [0] * (n)                                                                                            

            EulerPhi  $\gets$ [0, 1, 2, …, n]

            Mu  $\gets$ [1] * (n)

            Primes  $\gets$ [0] * (n)

            pos  $\gets$ 1

            for i = 2 to n do

if flag[i] = 0 then

                                    for j = i to n do

                                                flag[j]  $\gets$ 1

                                                EulerPhi[j]  $\gets$ ( EulerPhi[j] / i ) * ( i– 1 )

                                                if  j/i%i = 0 then

                                                            Mu[j]  $\gets$ 0

                                                else

                                                            Mu[j]  $\gets$ Mu[j] * (-1)

                                                end if

                                                j $\gets$ j + i

                                    end for

                                    Primes[pos]  $\gets$ i

                                    pos  $\gets$ pos +1

                           end if

          end for

          return Primes

          return EulerPhi

          return Mu

end function

This should be the end of this discussion. But I noted that it is as easy to initiate an array full of 1's (or any other value for that matter) as it is to initiate with 0. For example, when I tried to program the above algorithm in C, I first had to run the loop for the entire value of $n$, just to initialize the EulerPhi and Mobius function values. I therefore made a small optimization to the above algorithm which you can see here.

UPDATE (17/Nov/2019): Yet another interesting, and possibly faster, sieve is described in the Codeforces blog: Math note - linear sieve.

Yours Aye
Me