Cod sursa(job #2647555)

Utilizator MarcGrecMarc Grec MarcGrec Data 5 septembrie 2020 11:40:31
Problema Suma divizorilor Scor 80
Compilator cpp-64 Status done
Runda Arhiva de probleme Marime 2.66 kb
#define MOD 9901
#define MAX_P 7072

#include <fstream>
#include <vector>
#include <bitset>
#include <utility>
#include <cmath>
using namespace std;

ifstream fin("sumdiv.in");
ofstream fout("sumdiv.out");

int64_t a, b;
vector<int64_t> P;
vector<pair<int64_t, int64_t>> F;

void Erathosthenes(vector<int64_t>& out, int64_t limit);

void Factorization(const vector<int64_t>& Primes, vector<pair<int64_t, int64_t>>& out, int64_t num);

int64_t XNMOD(int64_t x, int64_t n, int64_t mod);

void Euclid(int64_t a, int64_t b, int64_t* d, int64_t* x, int64_t* y);

int64_t ModularInverse(int64_t num, int64_t mod);

int main()
{	
	fin >> a >> b;
	
	Erathosthenes(P, MAX_P);
	
	Factorization(P, F, a);
	
	int64_t prod1 = 1;
	int64_t prod2 = 1;
	for (const auto& factor : F)
	{
		int64_t aux = XNMOD(factor.first, b * factor.second + 1, MOD) - 1;
		
		while (aux < 0)
		{
			aux += MOD;
		}
		
		prod1 = (prod1 * aux) % MOD;
		
		prod2 = ((factor.first  - 1) * prod2) % MOD;
	}
	
	int64_t inverseProd2 = ModularInverse(prod2, MOD);
	
	int64_t divSumMod = (prod1 * inverseProd2) % MOD;
	
	fout << divSumMod;
	
	fin.close();
	fout.close();
	return 0;
}

vector<bool> prime;
void Erathosthenes(vector<int64_t>& out, int64_t limit)
{
	prime = vector<bool>(limit + 1, true);
	prime[0] = prime[1] = false;
	
	for (int64_t i = 2; (i * i) <= limit; ++i)
	{
		if (prime[i])
		{
			for (int64_t j = 2; (i * j) <= limit; ++j)
			{
				prime[i * j] = false;
			}
		}
	}
	
	out.clear();
	for (int64_t i = 0; i <= limit; ++i)
	{
		if (prime[i])
		{
			out.push_back(i);
		}
	}
	out.shrink_to_fit();
}

void Factorization(const vector<int64_t>& Primes, vector<pair<int64_t, int64_t>>& out, int64_t num)
{
	out.clear();
	for (int64_t prime : Primes)
	{
		if (prime > num)
		{
			break;
		}
		
		int64_t exp = 0;
		while ((num % prime) == 0)
		{
			++exp;
			num /= prime;
		}
		
		if (exp > 0)
		{
			out.emplace_back(prime, exp);
		}
	}
	
	if (num > 1)
	{
		out.emplace_back(num, 1);
	}
	
	out.shrink_to_fit();
}

int64_t XNMOD(int64_t x, int64_t n, int64_t mod)
{
	if (n == 0)
	{
		return 1;
	}
	
	int64_t res = XNMOD(x, n / 2, mod);
	
	res = (res * res * ((n & 1) ? x : 1)) % mod;
	
	return res;
}

void Euclid(int64_t a, int64_t b, int64_t* d, int64_t* x, int64_t* y)
{
	if (b == 0)
	{
		*d = a;
		*x = 1;
		*y = 0;
	}
	else
	{
		int64_t x0, y0;
		
		Euclid(b, a % b, d, &x0, &y0);
		
		*x = y0;
		*y = x0 - (a / b) * y0;
	}
}

int64_t ModularInverse(int64_t num, int64_t mod)
{
	int64_t x, y, d;
	
	Euclid(num, mod, &d, &x, &y);
	
	while (x < 0)
	{
		x += mod;
	}
	
	return x;
}