Cod sursa(job #997519)

Utilizator stefan.petreaStefan Petrea stefan.petrea Data 14 septembrie 2013 13:22:43
Problema Asmax Scor 0
Compilator cpp Status done
Runda Arhiva de probleme Marime 2.8 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
 *
 *
 *      r
 * A = Prod p_i^{a_i}
 *      i=1
 *
 *
 * 
 * theta_1(A^B) =  
 *
 *
 *     r   p_i^{a_i*B+1} - 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
// http://en.wikipedia.org/wiki/Modular_exponentiation
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;
};

/* There are probably many ways to compute the modular inverse. One is this:
 *
 *  if (a,m)=1 then
 *  a^-1 = a^phi(m)-1 (mod m)
 *
 *  http://mathworld.wolfram.com/ModularInverse.html
 *
 *  The other is by just using the extended euclidean algorithm
 *
 *  http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
 *
 *  Since I'm also computing a gcd, this is probably completely sub-optimal.
 *  
 */

inline long invmod(long a) {
    long d = gcd(a,mod);

    if(d==1) {
        return expmod(a,PHI-1);
    } else {
        //cout << "second\n"<<endl;
        return (expmod(d,PHI-1) * expmod(a/d,PHI-1)) % mod ;
    };
};

inline long compute_factor(long B,long a, long p) {
    /* compute one fraction from the big product */
    long v = invmod(p-1);
    a %= mod-1;
    long w = ((a * B) + 1) % mod;
    long u = expmod(p,w);
    //cout<<"a="<<a<<endl;
    //cout<<"B="<<B<<endl;
    //cout<<"w="<<w<<endl;
    //cout<<"u="<<u<<endl;
    u = u - 1;
    if(u < 0)
        u += mod;

    return  ((u % mod) * (v % mod)) % mod;
};


// powers of primes in A
int main() {
    ifstream I("sumdiv.in");
    ofstream O("sumdiv.out");
    long A,B;
    long p,a;
    long S = 1;
    I >> A >> B;

    a = 0;
    // Br = B % (mod-1)
    // it will be used to compute some powers % mod
    B %= mod-1;
    while(!(A & 1)) {
        A >>= 1;
        a++; 
    };
    S = ( S * compute_factor(B,a,2) ) % mod;

    for(p=3;p<=A;p+=2) {
        a = 0;

        //cout<<"p="<<p<<endl;
        if(A % p)
            return 1;

        while(A % p == 0) {
            A /= p ; 
            a++ ; 
        };

        S = (S * compute_factor(B,a,p)) % mod;
    };

    O << S;
    return 0;
};