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<j\le r}2\left\lfloor\frac{n}{p_ip_j}\right\rfloor

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

int answer=0;
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
			answer+=(ones*(m/p));
		else//there are an even number of 1s
			answer-=(ones*(m/p));
	}
}

This isn’t totally legal code because I had to paste it into HTML so the less than and greater than signs disappear and I don’t want to waste the time to fix it (even though I know how to).
For those of you that don’t know the bitset is just a string of 1s and 0s. When I did PIE=cnt I basically set it so that the bit string was equal to the base 2 representation of cnt. Using this method we can get all possible combinations of the prime factors.
The code can easily be converted to Java since I think it too has something similar to a bitset. If your language of choice doesn’t have a bit set, no need to worry. You can just get the base 2 representation on your own (it’s not that hard, trust me).
WARNING: This approach will only be feasible if each number doesn’t have too many prime factors. Not the case with the problem on SPOJ since each number can be decomposed into the product of at most 23 primes, and at most 9 distinct primes. So I’m not sure if the time limit will allow you use this method, I haven’t tried, but please let me know if you do.

Advertisements

One Response to Inclusion Exclusion

  1. hol says:

    I think you should also consider the cases where a prime factor is repeated, like values of 4, 8, 9. Your algorithm gives 35 for m=10 n=10, instead of 30. The extra prime counts are for (4,4),(4,8),(8,4),(8,8),(9,9), 5 in total. Is there some easy way to correct it? I coulnd’t figure out. Thanks.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: