Cod sursa(job #619782)

Utilizator claudiumihailClaudiu Mihail claudiumihail Data 16 octombrie 2011 04:50:07
Problema Principiul includerii si excluderii Scor 100
Compilator cpp Status done
Runda Arhiva educationala Marime 2.28 kb
#include <fstream>
#include <iostream>
#include <bitset>
#include <vector>
#include <cmath>

#define MAX_NUMS 1001001

using namespace std;

typedef unsigned long uint32;
typedef long long int64;
typedef unsigned long long uint64;

template<unsigned long numbits>
struct PrimeSieve
{
	uint64 num;
	bitset<numbits> bits;

	void GeneratePrimes(const uint64 n)
	{
		num=1;
		for (uint32 i=1; ((i<<1)+1)*((i<<1)+1)<n; ++i)
		{
			if (!bits[i])
			{
				for (uint64 j=i; ((i<<1)+1)*((j<<1)+1)<n; ++j)
				{
					const long long index=((i<<1)+1)*((j<<1)+1);
					bits[index>>1]=1;
				}
			}
		}
	}
	
	const bitset<numbits>& GetPrimes() const
	{
		return bits;
	}
};

uint32 divisors[MAX_NUMS];
PrimeSieve<16000002> primeSieve;

uint64 InclusionExclusion
(
	const uint64 A,
	const uint32 *divisors,
	const int numDivisors
)
{
	int idx=1;
	int64 sum=0;
	
	while (idx < (1<<numDivisors))
	{
		int aux = idx;
		int nr = -1;
		int64 prod = 1;

		while (aux)
		{
			const int index = __builtin_ffs(aux);
			prod *= divisors[index-1];
			nr++;
			aux = aux & (~(1<<(index-1)));
		}

		idx++;
		
		if (nr%2)
			nr = -1;
		else
			nr = 1;
		
		sum = sum + 1LL * nr * (uint64)(A/prod);
	}
	
	return sum;
}

int main()
{
	int m;
	uint64 a=1,b=1;
	vector<uint32> primes;
	fstream fin("pinex.in", fstream::in);
	fstream fout("pinex.out", fstream::out);

	primeSieve.GeneratePrimes(MAX_NUMS);

	primes.push_back(2);
	for(int i=1; (i<<1)+1<MAX_NUMS; ++i)
	{
		if(! (primeSieve.GetPrimes())[i] )
		{
			primes.push_back((i<<1)+1);
		}
	}
	
	fin >> m;

	for (int i=0; i<m; ++i)
	{
		fin >> a >> b;

		uint32 k=0, numDivs=0;
		uint32 sq = (int)ceil(sqrt(b));
		
		while (k<primes.size() && b>1)
		{
			if ((b%primes[k]) == 0)
			{
				divisors[numDivs] = primes[k];
				numDivs++;
				
				while ((b%primes[k]) == 0)
					b /= primes[k];
			}
			k++;
			
			// if we're over the sqrt and b is still > 1 it means there's one other prime that
			// can devide it...add it as the last prime divisor and exit
			if (primes[k] > sq && b>1)
			{
				divisors[numDivs] = b;
				numDivs++;
				break;
			}
		}

		fout << a - InclusionExclusion(a, divisors, numDivs) << "\n";
	}
	
	fin.close();
	fout.close();
	return 0;
}