Pagini recente » Cod sursa (job #3290070) | Cod sursa (job #1352991) | Cod sursa (job #997048)
Cod sursa(job #997048)
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
using namespace std;
#define min(x,y) (((x)<(y))?(x):(y))
#define max(x,y) (((x)>=(y))?(x):(y))
#define mod 9901
// euler's totient for 9901
#define PHI 9900
#define pL 10000
// first pL primes
/* the problem actually asks to compute theta_1(A^B) % 9901.
* http://en.wikipedia.org/wiki/Divisor_function
*
*
* theta_1(A^B) =
*
*
* r p_i^{a_i*B} - 1
* = Prod ---------------
* i=1 p_i - 1
*
* because 9901 is prime now, we can compute p_i^{a_i*B} fast through fermat's little theorem
*
* we also need the inverse mod 9901 of p_i - 1 which we can also compute fast through extended
* euclidean algorithm
*
*/
// compute a^b % mod
inline long long expmod(long long a,long long b) {
long long r = 1;
while(b > 0) {
if(b & 1)
r = (r * a) % mod;
b >>= 1;
a = (a*a) % mod;
};
return r;
};
inline long long gcd(long long a, long long b) {
long long r, i;
while(b!=0){
r = a % b;
a = b;
b = r;
}
return a;
};
/*
* if (a,m)=1 then
* a^-1 = a^phi(m)-1 (mod m)
*
* http://mathworld.wolfram.com/ModularInverse.html
*/
inline long long invmod(long long a) {
long long d = gcd(a,mod);
if(d==1) {
return expmod(a,PHI-1);
} else {
return (expmod(d,PHI-1) * expmod(a/d,PHI-1)) % mod ;
};
};
// powers of primes in A
int main() {
ifstream I("sumdiv.in");
ofstream O("sumdiv.out");
long long A,B;
long long i,j;
long long S = 1;
I >> A >> B;
// Br = B % (mod-1)
// it will be used to compute some powers % mod
long long Br = B % (mod-1);
for(i=2;i<=A;i++) {
long long p = i;
long long a = 0;
if(A % p)
continue;
while(A % p == 0) {
A /= p ;
a++ ;
};
/* compute one fraction from the big product */
long long u = 1;
long long v = invmod(p-1);
a %= mod-1;
long long w = ((a * Br) + 1) % mod;
u = expmod(p,w);
u = u - 1;
if(u < 0)
u += mod;
S = (S * ((u * v) % mod) ) % mod;
if(A == 1)
break;
};
O << S;
return 0;
};