Cod sursa(job #916550)
/**
suma divizorilor modulo 9901 lui A^B
1 <= A, B <= 50000000
SUMA = (1 + p1 + p1^2 + ... + p1^d1) * (1 + p2 + p2^2 + ... + p2^d2) * ... * (1 + pk + pk^2 + ... + pk^dk) =
= produs de progresii geometrice cu primul termen 1 si ratia p1, p2, ... ,pk =
(p1 ^ (d1+1) - 1)/(p1 - 1) * (p2 ^ (d2+1) - 1)/(p2 - 1) * ... * (pk ^ (dk+1) - 1)/(pk - 1)
descompun pe A in factori primi si la fiecare exponent inmultesc B si aplic formula de mai sus.
*/
#include <iostream>
#include <fstream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#define MOD 9901
using namespace std;
int A, B, sol, na, radA;
bool viz[8000];
int ciur[1500], nciur;
inline void CreareCiur()
{
radA = (int)sqrt(A*1.0);
ciur[++nciur] = 2;
int i, j, doii, radradA;
radradA = (int)sqrt(radA*1.0);
for (i=3; i<=radA; i+=2)
if (viz[i] == false)
{
ciur[++nciur] = i;
if (i <= radradA)
for (j=i*i, doii = i+i; j<=radA; j+=doii)
viz[j] = true;
}
}
inline int Putere (int x, int y)
{
/// x^y - 1
int put = 1;
while (y)
{
if (y&1)
{
put = (put * x%MOD)%MOD;
y--;
}
y >>= 1;
x = ((x%MOD)*(x%MOD))%MOD;
}
put--;
if (put<0)
put+=MOD;
return put;
}
inline int InversModular(int x)
{
/// x^(MOD-2)%MOD
int y, put;
put = 1;
y = MOD-2;
while (y)
{
if (y&1)
{
put = (put*(x%MOD))%MOD;
y--;
}
y >>= 1;
x = ((x%MOD)*(x%MOD))%MOD;
}
return put;
}
inline void Solve()
{
int i, div, put;
i = 1;
sol = 1;
while (i <= nciur && ciur[i]*ciur[i] <= A)
{
div = ciur[i];
i++;
if (A % div == 0)
{
put = 0;
while (A % div == 0)
{
put++;
A /= div;
}
put = put*B;
sol = (((sol * Putere(div, put+1))%MOD) * InversModular(div-1))%MOD;
}
}
if (A!=1)
{
if (A%MOD == 1)
{
/**
daca A%MOD == 1
sol = sol * (B+1) % MOD deoarece:
A^(B+1) - 1 = (A-1) * (A^B + A^(B-1) + ... + A + 1)
care se simplifica cu (A-1) de la numitor si obtinem
A^B + A^(B-1) + ... + A + 1 care are B+1 termeni si care
%MOD = B+1;
*/
sol = sol*((B+1)%MOD)%MOD;
}
else
sol = (((sol * Putere(A, B+1))%MOD) * InversModular(A-1))%MOD;
}
}
int main()
{
ifstream f("sumdiv.in");
f>>A>>B;
f.close();
CreareCiur();
Solve();
ofstream g("sumdiv.out");
g<<sol<<"\n";
g.close();
return 0;
}