## Inclusion Exclusion

I am sure that most of you have heard of the inclusion exclusion principle.  If not then I suggest that you either head over to Google or Wikipedia.  It is a really useful technique with many applications; even to this day I am still discovering new ways in which it can be used.

In this post I will outline a partial solution to the following problem located on SPOJ, PGCD.

The first thing you should notice is that the answer to (a,b) will be the same as the solution to (b,a).  This is because of the symmetry of the $gcd(\cdot,\cdot)$ operation, thus suppose that $a\le b$.

Suppose that there are $a$ columns in our gcd table.  We will find how many primes there are in each column.  Suppose that there are $m$ rows, and we are looking at column $n=p_1p_2\cdots p_r$, where $p_1,p_2,\cdots,p_r$ are distinct primes.  I’m not going to consider the case where a prime factor is repeated, though I have a feeling that the solution below should still work.

First we find how many multiples there are of each prime between $1\ \mbox{and}\ m$.  This is just

$\displaystyle\sum_{i=1}^r\left\lfloor\frac{n}{p_i}\right\rfloor$

But what about the multiples of say, $p_1p_2$?  Well they were counted twice, and their entry in the table is not going to be a prime (you can easily convince yourself of this fact).  So we just subtract them out twice, leaving our answer as

$\displaystyle\sum_{i=1}^r\left\lfloor\frac{n}{p_i}\right\rfloor-\sum_{1\le i

Now what about the divisors of $n$ with exactly 3 prime factors?  Their multiples have been removed too many times, thus we need to add them back.  If we continue this way we see that we will be using inclusion exclusion.

So to solve this we can use the following code.


#include bitset
#include vector

bitset PIE;
for(int i=2 to n;++i)
{
vector prime_factors=factor(i);
int r=prime_factors.size();
for(int cnt=1 to 2^r - 1;++cnt)
{
PIE=cnt;//express cnt in base 2 and store in PIE
int ones=PIE.count();//the number of 1s in our bit string
int p=1;
for(int j=0 to prime_factors.size() -1;++j)
if(PIE[j])
p*=prime_factors[j];
if(ones&1)//there are an odd number of 1s
else//there are an even number of 1s
}