Pagini recente » Rating Circiumaru David-Mihai (mihaidav) | Cod sursa (job #1383302) | Cod sursa (job #2045986) | Cod sursa (job #763463) | Cod sursa (job #997047)
Cod sursa(job #997047)
#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 expmod(long a,long b) {
long r = 1;
while(b > 0) {
if(b & 1)
r = (r * a) % mod;
b >>= 1;
a = (a*a) % mod;
};
return r;
};
inline long gcd(long a, long b) {
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 invmod(long a) {
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 A,B;
long i,j;
long S = 1;
I >> A >> B;
// Br = B % (mod-1)
// it will be used to compute some powers % mod
long Br = B % (mod-1);
for(i=2;i*i<=A;i++) {
long p = i;
long a = 0;
if(A % i)
continue;
while(A % p == 0) {
A /= p ;
a++ ;
};
/* compute one fraction from the big product */
long u = 1;
long v = invmod(i-1);
a %= mod-1;
long w = ((a * Br) + 1) % mod;
u = expmod(p,w);
u = u - 1;
if(u < 0)
u += mod;
S *= (u * v) % mod;
S %= mod;
if(A == 1)
break;
};
O << S;
return 0;
};