Cod sursa(job #996883)

Utilizator stefan.petreaStefan Petrea stefan.petrea Data 12 septembrie 2013 20:49:19
Problema Suma divizorilor Scor 10
Compilator cpp Status done
Runda Arhiva de probleme Marime 2.42 kb
#include <iostream>
#include <fstream>
#include <math.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; 
    int f = b&1;
    while(b > 1) {
        r = ((r % mod) * (r % mod)) % mod; 
        b/=2;
    };
    if(f)
        r = (r * b) % 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);
    } else {
        return (expmod(d,PHI) * expmod(a/d,PHI)) % 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<=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 */

        //cout<<"a["<<i<<"]="<<a<<endl;
        //cout<<"p["<<i<<"]="<<p<<endl;
        long u = 1;
        //cout<<"i="<<i<<endl;
        long v = invmod(i);
        a %= mod-1;
        long w = ((a * Br) + 1) % mod;
        //cout<<"Br="<<Br<<endl;
        //cout<<"w="<<w<<endl;
        for(j=1;j<=w;j++) {
            u *= p;
            u %= mod;
        };
        u = u - 1;
        if(u < 0)
            u += mod;
        //cout<<"u="<<u<<endl;
        //cout<<"v="<<v<<endl;
        S *= (u * v) % mod;
        S %= mod;

        if(A == 1)
            break;
    };

    O << S;
    return 0;
};