Cod sursa(job #482489)

Utilizator claudiumihailClaudiu Mihail claudiumihail Data 3 septembrie 2010 17:16:07
Problema Suma divizorilor Scor 80
Compilator cpp Status done
Runda Arhiva de probleme Marime 2.25 kb
#include<fstream>
#include<iostream>
#include<bitset>
#include<vector>
using namespace std;

typedef long long uint64;

const unsigned int MOD = 9901;

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

	void GeneratePrimes(const long long n)
	{
		num=1;
		for(int i=1; ((i<<1)+1)*((i<<1)+1)<n; ++i)
		{
			if(!bits[i])
			{
				for(long long 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;
	}
};

uint64 powl(uint64 n, int p)
{
	uint64 rez=1;
	while(p)
	{
		if(p&1)
		{
			rez=(rez*n)%MOD;
		}
		n=(n*n)%MOD;
		p>>=1;
	}
	return rez;
}

void Decompose(uint64& n, uint64 div, uint64& p, uint64& d)
{
	d=p=0;
	if(n % div == 0)
	{
		p=div;
		while(n % div == 0)
		{
			d++;
			n/=div;
		}
	}
}

int main()
{
	uint64 n,p,d,pow;
	fstream fin("sumdiv.in", fstream::in);
	fstream fout("sumdiv.out", fstream::out);
	vector<unsigned long> v;
	vector< pair<uint64, uint64> > factors;
	
	PrimeSieve<500007> primes;
	primes.GeneratePrimes(500005);
	v.push_back(2);
	for(int i=1; i<500005; ++i)
	{
		if(!(primes.GetPrimes())[i])
		{
			v.push_back((i<<1)+1);
		}
	}
	
	fin>>n>>pow;
	//cout<<n<<" ~ "<<pow<<endl;
	
	uint64 i=0;
	uint64 s=1,numdiv=1;
	while(n>1 && v[i]*v[i]<=n)
	{
		Decompose(n, v[i], p, d);
		if(p)
		{
			//cout<<p<<" - "<<d<<endl;
			factors.push_back(make_pair(p,d));
			/*numdiv*=(d+1);
			uint64 power=(powl(p, d+1)-1)%MOD;
			if(power<0)
				power+=MOD;
			uint64 divi=powl(p-1, MOD-2)%MOD;
			if(divi < 0)
				divi+=MOD;
			s=( ((s*power)%MOD)*divi)%MOD;*/
		}
		i++;
	}
	if(n>1)
	{
		//cout<<n<<" -- "<<1<<endl;
		//s=(s*(n+1))%MOD;
		factors.push_back(make_pair(n,1));
		numdiv*=2;
	}

	for(i=0; i<(int)factors.size(); ++i)
	{
		p=factors[i].first;
		d=pow*factors[i].second;
		//cout<<p<<" --- "<<d<<endl;
		uint64 power=(powl(p, d+1)-1)%MOD;
		if(power<0)
			power+=MOD;
		uint64 divi=powl(p-1, MOD-2)%MOD;
		if(divi < 0)
			divi+=MOD;
		s=( ((s*power)%MOD)*divi)%MOD;
	}
	fout<<s<<endl;
	//cout<<numdiv<<" "<<s<<"\n";
	
	fin.close();
	fout.close();
	return 0;
}