POJ1845 Sumdiv

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$ 的约数之和 mod 9901（1<=A,B<=5* $10^7$ ).
首先，我们需要熟练掌握因式分解。
$1+a+a^2+……+a^k$ = $(a^{k+1}-1)/(a-1)$
然后，根据算术基本定理：任何一个数都可以分为若干个（可以为1个）质数的幂之积。
由此可设：A= $p_1^{c_1} *p_2^{c_2} *……p_k^{c_k}$
所以A的约数之和为： $prod_{i=1}^k p_i^1+p_i^2+……p_i^{c_i}$ = $prod_{i=1}^k (p_i^{c_i+1}-1)/(p_i-1)$ mod 9901。
这个有点难理解，设A的约数x= $p_1^{d_1} *p_2^{d_2} *……p_k^{d_k}$ ，那么一定满足0<=di<=ci
假设 $p_1$ 的幂不确定，其他质数的幂之积y已经确定，那么约数和就为（ $p_1^0 +p_1^1 +……+p_1^{c_1}$ )*y.类似地，我们可以把y继续分解，就可以得到上述公式。（一定要弄懂，不懂可发问)
我们可设 $(p_i^{c_i}-1)=[(p_i-1)*x]*9901（x为正整数）+r（0<=r<9901*(p_i-1)）=[(p_i-1)*9901]*x+r$
那么 $(p_i^{c_i}-1)/(p_i-1)$ mod 9901= $r/(p_i-1)$ ，所以我们只要快速幂求r就行。
当然，还有一种做法，就是把 $(p_i^{c_i}-1)$ $mod$ $9901$ ,和 $(p_i-1)$ $mod$ $9901$ 求出，然后用乘法逆元的方法求出商（用到了费马小定理，因为 $a/b≡a * b^{9901-2}$ (mod 9901)（当 $gcd(b,9901)=1$ ))。当gcd(b,9901)!=1时，根据裴蜀定理（即把式子转化成同余方程），不能乘法逆元（得不出整数解），此时b≡1(mod 9901)（即 $p_i$ $mod$ $9901$ =1, $p_i^1+p_i^2+……p_i^{c_i}$ mod 9901= $c_i$ +1)。
以上我们求的是A的约数之和， $A^B$ 的约数之和只需把 $c_i$ 均乘以B就行。
下面贴两个代码：

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define ll long long
const ll mod=9901;
ll x,y,tot,s[33],b[33];
// x,y为问题中的a,b
//tot为质因数种数，s[i]表示第i种质数的个数，b[i]为第i种质数的值
void divide(ll k)
{
	for(ll i=2;i*i<=k;i++)
	{
		if(!(k%i))
		{
			b[++tot]=i;
			do
			{
				++s[tot];
				k/=i;
			}while(!(k%i));
		}
	}
	if(k!=1)b[++tot]=k,++s[tot];
}
void mult(ll &a,ll b,ll p)//a*b%p 两种方法 防溢出。注意出题人会千方百计出溢出的数据。 
{
	a%=p;b%=p;
	ll c=(long double)a*b/p;
	a=a*b-c*p;
	/*ll ans=0;a%=p;b%=p;
	while(b)
	{
		if(b)ans+=a,ans%=p;
		b>>=1;a<<=1;a%=p;
	}
	return ans;*/
}
ll power(ll a,ll b,ll c)//a^b%c
{
	ll ans=1%c;a%=c;
	while(b)
	{
		if(b&1)mult(ans,a,c);
		b>>=1;mult(a,a,c);
	}
	return ans;
}

int main()
{
	scanf("%lld%lld",&x,&y);
	if(x==1){puts("1");return 0;}
	else if(x==0){puts("0");return 0;}//特判 
	divide(x);
	ll ans=1;
	while(tot)
	{
		ll i=b[tot];
		ans*=((power(i,(y*s[tot--]+1),(i-1)*mod)-1+(i-1)*mod)%((i-1)*mod))/(i-1);
		ans%=mod;
	}
	printf("%lld
",ans);
	return 0;
}

//不知道为什么这样更快，可能是算模数更快吧。
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define ll long long
const ll mod=9901;
ll x,y,tot,s[33],b[33];
// x,y为问题中的a,b
//tot为质因数种数，s[i]表示第i种质数的个数，b[i]为第i种质数的值
void divide(ll k)
{
	for(ll i=2;i*i<=k;i++)
	{
		if(!(k%i))
		{
			b[++tot]=i;
			do
			{
				++s[tot];
				k/=i;
			}while(!(k%i));
		}
	}
	if(k!=1)b[++tot]=k,++s[tot];
}
void mult(ll &a,ll b,ll p)//a*b%p 两种方法 防溢出 
{
	a%=p;b%=p;
	ll c=(long double)a*b/p;
	a=a*b-c*p;
	/*ll ans=0;a%=p;b%=p;
	while(b)
	{
		if(b)ans+=a,ans%=p;
		b>>=1;a<<=1;a%=p;
	}
	return ans;*/
}
ll power(ll a,ll b,ll c)//a^b%c
{
	ll ans=1%c;a%=c;
	while(b)
	{
		if(b&1)mult(ans,a,c);
		b>>=1;mult(a,a,c);
	}
	return ans;
}

int main()
{
	scanf("%lld%lld",&x,&y);
	if(x==1){puts("1");return 0;}
	else if(x==0){puts("0");return 0;}//特判 
	divide(x);
	ll ans=1;
	while(tot>0)
	{
		ll i=b[tot];//取质数 
		if((i-1)%mod==0)ans*=s[tot]*y+1;
		else 
		{
			ll aa=(power(i,s[tot]*y+1,mod)-1+mod)%mod,bb=power(i-1,mod-2,mod);//求逆元 
			ans*=aa*bb;
		}
		ans%=mod;--tot;
	}
	printf("%lld
",ans);
	return 0;
}