【数论】Sumdiv(求X^Y因子和)

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

题意：
求X^Y的因子和，结果取模9901.

思路分析：
就是求因子和的模板题，运用到素数筛，合数分解，同余定理，快速乘法，快速幂，求因子和的公式。做这道题有几个坑点：

1. 需要运用快速乘法优化快速幂，不然会爆long long
2. 这道题不能用逆元做，虽然模9901是质数，但当a为9901的倍数时就无法互质。需要运用同余定理的除法通式。

下面直接放AC代码。

#include <iostream>
#include <cstring> 
#include <cmath>
#define ll long long
using namespace std; 

const int N = 200000;
ll prime[N];
int vis[N];
int t;

ll fast_mul(ll x,ll y,ll mod){	//快速乘法(x*y取模，防止溢出并且加快速度) 
    ll res=0;
    while(y){
        if(y&1) res=(res+x)%mod;
        x=(x+x)%mod;
        y /= 2;
    }
    return res;
}

ll pow_mod(ll x,ll y,ll mod){	//快速乘法优化快速幂 
    ll res=1;
    while(y){
        if(y&1) res=fast_mul(res,x,mod);
        x=fast_mul(x,x,mod);
        y /= 2;
    }
    return res;
}

void isPrime()	//素数筛 
{
    t=0;
    memset(vis,0,sizeof(vis));
    memset(prime,0,sizeof(prime));
    for(ll i=2;i<N;i++)
    {
        if(!vis[i])
        {
            prime[t++]=i;
            for(ll j=i*i;j<N;j+=i)
                vis[j]=1;
        }
    }
}

//factor[i][0]=pi,factor[i][1]=ai;
ll factor[100][2];
ll cnt;	//cnt表示质因数个数 

void splitPrime(ll x)	//分解质因数 
{
    cnt=0;
    ll t=x;
    for(int i=0;prime[i]<=t/prime[i];i++)
    {
        factor[cnt][1]=0;	 
        while(t%prime[i]==0)	//可以被质数整除 
        {
            factor[cnt][0]=prime[i];
            while(t%prime[i]==0)	//计算出除以质数的次数当作该质数的幂 
            {
                factor[cnt][1]++;
                t/=prime[i];
            }
            cnt++;
        }
    }
    if(t!=1)	//t!=1表示t为质数 
    {
        factor[cnt][0]=t;
        factor[cnt][1]=1;
        cnt++;
    }
}

ll sum_factor(ll x, ll y, ll c) //求因子总和。s(n) = (1+p^1+p^2+...p^e)*s(n') = (p^(e+1) - 1) / (p-1) * s(n')
{								  //令g(p, e) = (p^(e+1)-1) / (p-1)，则s(n) = g(p1, e1) * g(p2, e2) * ... * g(pk, ek)
	if(x == 0) return 0;
	splitPrime(x);
	ll sum = 1;
	for(int i = 0;i < cnt;i++)
	{
		ll b = factor[i][0]-1;
		ll m = b*c;
		ll a = pow_mod(factor[i][0],y*factor[i][1]+1,m)-1;
		ll tmp = a/b%c;
		sum *= tmp%c; 
		sum = (sum+c)%c; 
	}
	return sum;
} 

int main()
{
	isPrime();
	ll A,B;
	while(~scanf("%lld%lld",&A, &B))
	{
		printf("%lld\n",sum_factor(A,B,9901));
	}
	return 0;
}

补充：关于对tmp算法的解释