POJ 1845 Sumdiv（逆元的应用）

传送门

Sumdiv

Time Limit: 1000MS		Memory Limit: 30000K
Total Submissions: 19009		Accepted: 4773

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).

Source

Romania OI 2002

题目大意：

求 $A^B$ 的因子和对 $MOD$ 取模。

解题思路：

首先　对 $A$ 进行素因子分解：

$A = p_1^{a_1} *p_2^{a_2}*...*p_k^{a_k}$

素因子分解完之后 $A^B$ 就等于

$$p_1^{a_1*B} *p_2^{a_2*B}*...*p_k^{a_k*B}$$
然后他们的素因子的和就可以写作：
$$ans = (1+p_1+p_1^2+...+p_1^{a_1*B})*((1+p_2+p_2^2+...+p_2^{a_2*B}))*...*(1+p_k+p_k^2+...+p_k^{a_k*B})$$
每一个括号中的数是一个等比数列，然后等比数列求和得到：
$$ans = \frac{p_1^{a_1*B+1}}{p_1-1}*\frac{p_2^{a_2*B+1}}{p_2- 1}*...\frac{p_k^{a_k*B+1}}{p_k-1}$$
现在我们要求的就是 $ans \% MOD$ ，但是 $ans$ 带着分数，那么就要求逆元，可是这里有一个问

题是 逆元可能不存在，那么怎么办呢，我们就用到了一个公式：

$$\frac AB\%MOD =\frac{A\%(MOD*B)}B$$
根据这个公式就可以做了，首先要素数筛，然后素因子分解，快速幂，这里需要注意的是，当我们

进行快速幂的时候，中间有可能爆 $long$ $long$ 所以，在快速幂的过程中加上快速乘法，这样就可

以了。

$My$ $Code：$

/**
2016 - 08 - 08 下午
Author: ITAK

Motto:

今日的我要超越昨日的我，明日的我要胜过今日的我，
以创作出更好的代码为目标，不断地超越自己。
**/

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <vector>
#include <queue>
#include <algorithm>
#include <set>
using namespace std;
typedef long long LL;
typedef unsigned long long ULL;
const int INF = 1e9+5;
const int MAXN = 1e6+5;
const LL MOD = 9901;
const double eps = 1e-7;
const double PI = acos(-1.0);
using namespace std;
LL p[MAXN];
int k = 0;
bool prime[MAXN];
void isprime()
{
    k = 0;
    memset(prime, 0, sizeof(prime));
    prime[1] = 1;
    for(int i=2; i<MAXN; i++)
    {
        if(!prime[i])
        {
            p[k++] = i;
            for(int j=i+i; j<MAXN; j+=i)
                prime[j] = 1;
        }
    }
}
LL fac[MAXN/100], num[MAXN/100];
int cnt = 0;
void Dec(LL x)
{
    cnt = 0;
    memset(num, 0, sizeof(num));
    for(int i=0; p[i]*p[i]<=x; i++)
    {
        if(x%p[i] == 0)
        {
            fac[cnt] = p[i];
            while(x%p[i] == 0)
            {
                num[cnt]++;
                x /= p[i];
            }
            cnt++;
        }
    }
    if(x > 1)
    {
        fac[cnt] = x;
        num[cnt++] = 1;
    }
}
LL multi(LL a, LL b, LL c)
{
    LL ans = 0;
    while(b)
    {
        if(b & 1)
            ans = (ans+a)%c;
        b>>=1;
        a = (a+a)%c;
    }
    return ans;
}
LL quick_mod(LL a, LL b, LL c)
{
    LL ans = 1;
    while(b)
    {
        if(b & 1)
            ans = multi(ans, a, c);///可能爆 long long
        b>>=1;
        a = multi(a, a, c);
    }
    return ans;
}

int main()
{
    isprime();
    LL A, B;
    while(~scanf("%lld%lld",&A,&B))
    {
        Dec(A);
        LL ans = 1, x, y, tmp;
        for(int i=0; i<cnt; i++)
        {
            LL MM = MOD*(fac[i]-1);///MM * MM 容易爆 long long
            tmp = quick_mod(fac[i], num[i]*B+1, MM);
            tmp--;
            tmp = (tmp%MM+MM)%MM;
            tmp /= (fac[i]-1);
            ans = (ans*tmp)%MOD;
        }
        ans = (ans%MOD+MOD)%MOD;
        printf("%lld\n",ans);
    }
    return 0;
}