POJ-1845 Sumdiv (快速乘)(乘法逆元)(唯一分解定理)

題幹：

給定a、b，求 a b a^b ab的所有因子和，結果模9901.

0 <= A,B <= 50000000

思路：

根據唯一分解定理：a= P 1 e 1 P1^{e1} P1e1 * P 2 e 2 P2^{e2} P2e2 * P 3 e 3 P3^{e3} P3e3 … P n e n Pn^{en} Pnen（Pi為素數）

是以a的因子和可以寫成：ans=(1+ P 1 1 P1^{1} P11+ P 1 2 P1^{2} P12… P 1 e 1 P1^{e1} P1e1) * (1+ P 2 1 P2^{1} P21+ P 2 2 P2^{2} P22… P 2 e 2 P2^{e2} P2e2) …*(1+ P n 1 Pn^{1} Pn1+ P n 2 Pn^{2} Pn2… P n e n Pn^{en} Pnen)

是以 a b a^b ab的因子和可以寫成：ans=(1+ P 1 1 ∗ b P1^{1*b} P11∗b+ P 1 2 ∗ b P1^{2*b} P12∗b… P 1 e 1 ∗ b P1^{e1*b} P1e1∗b) * (1+ P 2 1 ∗ b P2^{1*b} P21∗b+ P 2 2 ∗ b P2^{2*b} P22∗b… P 2 e 2 ∗ b P2^{e2*b} P2e2∗b) …*(1+ P n 1 ∗ b Pn^{1*b} Pn1∗b+ P n 2 ∗ b Pn^{2*b} Pn2∗b… P n e n ∗ b Pn^{en*b} Pnen∗b)

觀察(1+ P 1 1 ∗ b P1^{1*b} P11∗b+ P 1 2 ∗ b P1^{2*b} P12∗b… P 1 e 1 ∗ b P1^{e1*b} P1e1∗b) 可得這是一個等比數列的和，可以寫成

（ P 1 e 1 ∗ b + 1 P1^{e1*b+1} P1e1∗b+1-1）/（P1-1）%9901.

因為除法不能取餘，是以我們對（P1-1）取乘法逆元，變成：

（ P 1 e 1 ∗ b + 1 P1^{e1*b+1} P1e1∗b+1-1）-1%(9901*(P1-1))/(P1-1)

因為牽扯到減法是以變成：

（ P 1 e 1 ∗ b + 1 P1^{e1*b+1} P1e1∗b+1-1)+(9901*(P1-1))-1%(9901*(P1-1))/(P1-1)

然後對 P 1 e 1 ∗ b + 1 P1^{e1*b+1} P1e1∗b+1做快速幂，為了防止乘的中途爆資料，使用快速乘(即将一個乘數按二進制乘然後相加)。

#include <cstdio>  
#include <cstring>
#include <iostream>    
#include <algorithm>
#include <cmath>
#include <map>
using namespace std;
typedef long long ll;
const ll mod=9901;
const int mx=50005;
int f[mx],p[mx],num;
void f1()
{
	for(int i=2;i<mx;i++){
		if(!f[i]) p[num++]=i;
		for(int j=0;j<num&&i*p[j]<mx;j++){
			f[i*p[j]]=1;
			if(!(i%p[j]))  break;
		}
	}
}
ll mul(ll a,ll b,ll c)
{
	a%=c;
	b%=c;
	ll ans=0;
	while(b)
	{
		if(b&1)
			ans=(ans+a)%c;
		a=(a+a)%c;
		b>>=1;
	}
	return ans%c;
}
ll qc(ll a,ll b,ll c)
{
	ll ans=1;
	while(b)
	{
		if(b&1)
			ans=mul(ans,a,c);
		a=mul(a,a,c);
		b>>=1;
	}
	return ans%c;
}
int main()
{
	ll a,b;
	f1();
	scanf("%lld %lld",&a,&b);
	ll ans=1;
	for(int i=0;p[i]*p[i]<=a;i++){
		int sum=0;
		if(a%p[i]==0)
		{
			while(a%p[i]==0)
			{
				sum++;
				a/=p[i];
			}
			ll m=(p[i]-1)*mod;
			ans*=(qc(p[i],sum*b+1,m)+m-1)/(p[i]-1);
			ans%=mod;
		}
		//printf("%lld\n",ans);
	}
	if(a>1)
	{
		ll m=(a-1)*mod;
		ans*=(qc(a,b+1,m)+m-1)/(a-1);
		ans%=mod;
	}
	printf("%lld\n",ans);
	return 0;
}