POJ 1845 逆元最清真写法

Description

Consider two natural numbers A and B. Let S be the sum of all natural divisors of A^B. Determine S modulo 9901 (the rest of the division of S by 9901).

Input

The only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.

Output

The only line of the output will contain S modulo 9901.

Sample Input

2 3

Sample Output

15

Hint

2^3 = 8.
The natural divisors of 8 are: 1,2,4,8. Their sum is 15.
15 modulo 9901 is 15 (that should be output).

---

就是求$A^B$的所有因数之和

---

对于一个数$A$
唯一分解为$A=a_1^{k_1}*a_2^{k_2}*a_3^{k_3}*a_4^{k_4}*...*a_n^{k_n}$
则约数和为
$Sum=(a_1^{0}+a_1^{1}+a_1^{2}+...+a_1^{k_1})*(a_2^{0}+a_2^{1}+a_2^{2}+...+a_2^{k_2})*...*(a_n^{0}+a_n^{1}+a_n^{2}+...+a_n^{k_n})$

---

注意以下两个点
①：如果分解的某个质数-1之后mod 9901 为0，则这个东西的贡献等同于它的次数+1
②：要（%mod+mod）%mod

#include<cstdio>
#include<iostream>
#include<cstring>
using namespace std;
const int mod = 9901;
long long p[1000000], cnt[1000000];
long long fast_pow( long long x, long long k ) {
    long long te = 1, mp = x;
    while( k ) {
        if( k & 1 ) te = te * mp % mod;
        mp = mp * mp % mod;
        k >>= 1;
    }
    return te % mod;
}
void exgcd(long long a, long long b, long long& d, long long& x, long long& y ) {
    if ( b == 0 ){ d = a; x = 1; y = 0; }
    else { exgcd ( b, a%b, d, y, x); y -= x * ( a / b ); }
}
long long inverse ( long long a, long long b ) {
    long long d, x, y;
    exgcd( a, b, d, x, y);
    return ( x % mod + mod ) % mod;
}
long long sum( long long p, long long k ){
    long long tmp=inverse( p-1, mod );
    return ( p * ( fast_pow ( p, k ) - 1 ) * tmp +1 ) % mod;
}
int main() {
    long long A, B;
    cin >> A >> B;
    if ( A == 0 ) {
        printf( "0\n" ); return 0;
    }
    int k = 0;
    for( register long long i = 2; i * i <= A; ) {
        if( A % i == 0 ){
            p[++k] = i;
            cnt[k] = 0;
            while ( ! ( A % i ) ) cnt[k]++ , A /= i;
        }
        if( i == 2) i++;
        else        i+=2;
    }
    if( A != 1){
        p[++k] = A;
        cnt[k] = 1;
    }
    long long ans = 1;
    for( register long long i = 1; i <= k; i++){
        if( ( p[i] - 1 ) % mod == 0 ) ans = ( ans * ( cnt[i] * B % mod + 1 ) % mod + mod ) % mod;
        else ans = ( ans * sum ( p[i], cnt[i] * B ) % mod + mod )% mod;
    }
    cout << ans << endl;
    return 0;
}