Sunday, November 29, 2015

C-code for the PrimePi function

This code uses the algorithm described here. The same algorithm has been used for SumOfPrimes function here.


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

#define M 1000000
#define Mod 1000000000

int Lower[M+1], Higher[M+1];
char Used[M+1];

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 PrimePi(long long n)
{
    long long v = sqrt(n), p, temp, q, j, end, i, d, t;
    Higher[1] = n % Mod - 1;
    for(p = 2; p <= v; p++) {
        Lower[p] = p - 1;
        Higher[p] = (n/p) % Mod - 1;
    }
    for(p = 2; p <= v; p++) {
        if(Lower[p] == Lower[p-1])
            continue;
        temp = Lower[p-1];
        q = p * p;
        Higher[1] -= Higher[p] - temp;
        if(Higher[1] < 0)
            Higher[1] += Mod;
        j = 1 + (p & 1);
        end = (v <= n/q) ? v : n/q;
        for(i = p + j; i <= 1 + end; i += j) {
            if(Used[i])
                continue;
            d = i * p;
            if(d <= v)
                Higher[i] -= Higher[d] - temp;
            else {
                t = n/d;
                Higher[i] -= Lower[t] - temp;
            }
            if(Higher[i] < 0)
                Higher[i] += Mod;
        }
        if(q <= v)
            for(i = q; i <= end; i += p*j)
                Used[i] = 1;
        for(i = v; i >= q; i--) {
            t = i/p;
            Lower[i] -= Lower[t] - temp;
        }
    }
    return Higher[1] % Mod;
}

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

Well, I conclude the discussion on the PrimePi functions with this post, I'll think of something to continue blogging which may take a while.



Until then
Yours Aye
Me

C-code for the SumOfPrimes function


We saw the algorithm for the generalized Prime counting function in the previous post. Here is the C code for the same, with some slight optimizations.


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

#define M 1000000
#define Mod 1000000000

long long Lower[M+1], Higher[M+1];
char Used[M+1];

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 SumPrimePi(long long n)
{
    long long v = sqrt(n), p, temp, q, j, end, i, d, t;
    Higher[1] = SumNat(n) - 1;
    for(p = 2; p <= v; p++) {
        Lower[p] = SumNat(p) - 1;
        Higher[p] = SumNat(n/p) - 1;
    }
    for(p = 2; p <= v; p++) {
        if(Lower[p] == Lower[p-1])
            continue;
        temp = Lower[p-1];
        q = p * p;
        Higher[1] -= (Higher[p] - temp) % Mod * p % Mod;
        if(Higher[1] < 0)
            Higher[1] += Mod;
        j = 1 + (p & 1);
        end = (v <= n/q) ? v : n/q;
        for(i = p + j; i <= 1 + end; i += j) {
            if(Used[i])
                continue;
            d = i * p;
            if(d <= v)
                Higher[i] -= (Higher[d] - temp) % Mod * p % Mod;
            else {
                t = n/d;
                Higher[i] -= (Lower[t] - temp) % Mod * p % Mod;
            }
            if(Higher[i] < 0)
                Higher[i] += Mod;
        }
        if(q <= v)
            for(i = q; i <= end; i += p*j)
                Used[i] = 1;
        for(i = v; i >= q; i--) {
            t = i/p;
            Lower[i] -= (Lower[t] - temp) % Mod * p % Mod;
            if(Lower[i] < 0)
                Lower[i] += Mod;
        }
    }
    return Higher[1] % Mod;
}

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


The code computes the value of $\pi_1(10^{12})$ in about 6 seconds which is way, way better than Mathematica's computation time. If you are looking for the C-code for $\pi(x)$, you can find it in this post.

Saturday, November 28, 2015

Algorithm for the Prime counting function


In Counting Primes - Prime Counting function - Part III, we saw the function used to calculate the sum of all primes less than a given value and it's Mathematica code. Well, the Mathematica code is good enough but we can do better with C. But if we try to convert the code given into C, we run into a problem. In the given code, we have relied heavily upon Mathematica's in-built functions. But we can sort the problem easily by using the some of the properties $\varphi_k(x,\pi(\sqrt{x}))$. First note that,

Result 1 : $\pi_k(x)=\varphi_k(x,\pi(\sqrt{x}))=\varphi_k(x,\sqrt{x})=\varphi_k(x,x-2)=\varphi_k(x,x-1)$, since $\pi(x) \leq x$.

Secondly, $\pi_k(x)$ increases only at prime values of $x$. This can be used to test the primality of $x$ in our C-code. Or casting it in different terms, we can say

Result 2 : $x$ is prime, if and only if $\pi_k(x) > \pi_k(x-1)$ (or) $\varphi_k(x,x-2) > \varphi_k(x-1,x-2)$.



The Algorithm to find $\pi_k(n)$ will now proceed as follows.

Input : Get the value of $n$.

Step 1: Create a table with 2 rows. Mark the first row as $m$ and fill it as $1,2,\cdots,\lfloor\sqrt{n}\rfloor, \lfloor n/\lfloor\sqrt{n}\rfloor\rfloor,\cdots,\lfloor n/2\rfloor,\lfloor n/1 \rfloor$. Mark the second row as $\varphi_k(m)$ and the fill the corresponding $\varphi_k(m,0)$ values.
Step 2 : Initiate $r=2$.
Step 3 : Check whether $r$ prime or not. To do this, compare the corresponding $\varphi$ entries of $r$ and $r-1$ (making use of Result 2).
Step 4 : If they are same, skip to Step 6. Else, conclude $r$ to be a prime.
Step 5 : For each $m\geq r^2$ , starting from the last entry and moving left, update the values of $\varphi_1$ using the recurrence, $\varphi_k(m)=\varphi_k(m)-r^k*(\varphi_k(\lfloor m/r \rfloor)-\varphi_k(r-1))$
Step 6 : Increment $r$.
Step 7: If $r$ equals $\lfloor\sqrt{n}\rfloor$, output the rightmost entry of the second row and Stop. Else, goto Step 3.



For example, the tables while computing $\pi_1(20)$ are given below.

TABLE 1 :: $\varphi_1(m,0)$
$
\begin{matrix}
m&1&2&3&4&5(=\lfloor20/4\rfloor)&6(=\lfloor20/3\rfloor)&10(=\lfloor20/2\rfloor)&20(=\lfloor20/1\rfloor)\\
\varphi_1(m)&0&2&5&9&14&20&54&209\\
\end{matrix}
$

TABLE 2 :: $\varphi_1(m,1)$
$
\begin{matrix}
m&1&2&3&4&5&6&10&20\\
\varphi_1(m)&0&2&5&5&10&10&26&101\\
\end{matrix}
$

TABLE 3 :: $\varphi_1(m,2)$
$
\begin{matrix}
m&1&2&3&4&5&6&10&20\\
\varphi_1(m)&0&2&5&5&10&10&17&77\\
\end{matrix}
$

Note a couple of things here. The table is modified $\pi(\lfloor\sqrt{n}\rfloor)$ times after the initialization. As a byproduct of the algorithm, we get the SumOfPrimes for all the entries in the first row. Also, we generate the list of primes $\leq \lfloor\sqrt{n}\rfloor$.

After I wrote this, I found that no text can really equal a cute demonstration in terms of clarity. So I found this, which will visualizes the algorithm's execution. Hope that clarifies things. I'll post the C-code for this algorithm in the next post.


Until then
Yours Aye
Me

Friday, November 27, 2015

Pólya's Enumeration Theorem with Applications - Part II

This post continues what we discussed here.

Pólya's Enumeration theorem allows us to do much more. In all of the examples above, we did not care about the number of balls we had for each color. In fact, that information was not needed to answer the questions we asked ourselves. Also, in both the necklace and the urn problems, we allowed only one ball at a given position. What if there can be more than one ball at one particular position? Consider the following question.

Example 5: 8 indistinguishable Urns arranged in a circle. What is the number of ways to distribute $n$ blue balls into these urns discounting rotational symmetry?

To make the question clear, number the urns from 1 to 8 in clockwise direction. Let's say we have 16 blue balls. Then according to the question, the distributions {1,2,3,4,5,1,0,0} and {5,1,0,0,1,2,3,4} are identical because one can be transformed to other by rotations. To solve this problem with Pólya's Enumeration Theorem we first find the cycle index of Cyclic group of degree 8.

$Z_{C_8}(x_1,x_2,...,x_8)=\frac{1}{8}(x_1^8+x_2^4+2x_4^2+4x_8)$

Pólya's Enumeration theorem says the answer we are looking for is given by

$\frac{1}{8}((X_1^8)_n+(X_2^4)_n+2(X_4^2)_n+4(X_8)_n)$

where $X_a^b=\frac{1}{(1-z^a)^b}$ and $(F(z))_n$ is the coefficient of $z^n$ in the power series of $F(z)$. With Mathematica's help we find the answer to be $30667$ when $n=16$.

The same procedure can be applied for all the permutations groups and each will give the required answer according to their nature. One thing of interest here is that when there are $n$ balls of the same color, $Z_{S_r}$ counts the number of partitions of $n$ into $r$ parts ($0$ included) and $Z_{-S_r}$ counts the number of partitions of $n$ into $r$ distinct parts ($0$ included). This interpretation has an interesting application in the multiplicative partition of $n$ which will be final piece of our discussion.

Now, the previous problem had balls of the same color. What if there are $n_i$ balls of color $c_i$. This is the second main point that I wished to convey in this post. Let's look at the next problem.

Example 6: 4 indistinguishable Urns are arranged in a circle. What is the number of ways to distribute $m$ different colored balls, with $n_i$ balls of color $c_i$, into these urns discounting rotational symmetry?

The answer to this question proceeds in pretty much the same way as the previous answer. We find the cycle index of Cyclic group of degree 4. The important distinction comes in the definition of $(F(z))_n$. Write $n=(n_1,n_2,..,n_m)$ and define,

$(F(z))_n=\Pi_{j=1}^m(F(z))_{n_j}$ where $F(z)$ will be of the form $\Pi_{i=1}^rX_i^\alpha$.

In case of 4 Urns with 1 white balls and 1 black balls discounting rotational symmetry, we use

$Z_{C_4}(x_1,x_2,x_3,x_4)=\frac{1}{4}(x_1^4+x_2^2+2x_4)$

and the final answer will be $\frac{1}{4}((X_1^4)_{(1,1)}+(X_2^2)_{(1,1)}+2(X_4)_{(1,1)}) = 4$

because {bg,$0$,$0$,$0$}, {b,$0$,g,$0$}, {b,g,$0$,$0$} and {g,b,$0$,$0$} are the possible arrangements.

Example 7: What is the number of ways to distribute 2 white balls and 1 black balls into 3 indistinguishable Urns without taking order into account?

This question can be interpreted in at least two (equivalent) ways. First, it is the number of ways of partitioning the 2-partite number (2,1) into three 2-partite numbers. Second, it is the number of multiplicative partitions of $12(=2^2\times3^1)$. That's easy. $12 = 1\times 2\times 6 = 1\times 3\times 4 = 2\times 2\times 3 = 1\times 1\times 12$. Let's check whether Pólya's Enumeration Theorem gives us the same answer.

$Z_{S_3}=\frac{1}{6}(x_1^3+3x_1x_2+2x_3)$

Following the same procedure as above we see the answer is $\frac{1}{6}((X_1^3)_{(2,1)}+3(X_1X_2)_{(2,1)}+2(X_3)_{(2,1)})$.  Let's just evaluate one of these terms, preferably $(X_1X_2)_{(2,1)}$.

$(X_1X_2)_{(2,1)}=[z^2]\frac{1}{(1-z)(1-z^2)}\times[z^1]\frac{1}{(1-z)(1-z^2)}=2\times1=2$.

The final answer is $\frac{1}{6}(6\times3+3\times2\times1+2\times0\times0)=4$ which we saw already. In terms of balls and Urns the distribution would look like {$0$,w,wb}, {$0$,b,ww}, {w,w,b} and {$0$,$0$,wwb}.

Example 8: What is the answer to the same question if we want no two urns to have the same pattern of balls?

In this case it asks for the multiplicative partition of 12 with distinct numbers in which case the last two multiplicative partitions listed above will be eliminated. In terms of Urns and balls, {w,w,b} will be eliminated because two urns have one white ball each and {$0$,$0$,wwb} will be eliminated because two urns are empty. To get the answer via Pólya's Enumeration theorem, we'll be using $Z_{-S_3}$ instead of $Z_{S_3}$.

$Z_{-S_3}=\frac{1}{6}(x_1^3-3x_1x_2+2x_3)$

Hence the answer will be  $\frac{1}{6}((X_1^3)_{(2,1)}-3(X_1X_2)_{(2,1)}+2(X_3)_{(2,1)})=2$.

Well, that's it from my end. I'm quite impressed with the way I progressed but it's still upto a first time reader to really judge how good I was in explaining what I intended to explain. Any suggestions and corrections are welcome. It feels extremely good to have come up with this post and I enjoyed it. Hope you'll do too. C ya later.

Until then
Yours Aye,
Me

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

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

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