Cod sursa(job #997048)

Utilizator stefan.petreaStefan Petrea stefan.petrea Data 13 septembrie 2013 08:59:46
Problema Suma divizorilor Scor 40
Compilator cpp Status done
Runda Arhiva de probleme Marime 2.24 kb
#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;
};