poj1845 Sumdiv(算数基本定理,质因数分解)

poj

(1<=a,b<=5*10^7))

((1+p_1+p_1^2+...+p1^{b*c_1})*...*(1+p_n+p_n^2+...+pn^{b*c_n}))

(1+p_1+p_1^2+...+p1^{b*c_1}=(p_1^{b*c_1+1}-1)/(p_1-1))

(p_1^{b*c_1+1}-1).

(1+p_1+p_1^2+...+p1^{b*c_1}=1+1+1^2+..+1^{b*c_1}=b*c_1+1),也是可以直接计算的.

const int mod=9901;
int a,b,m,ans=1;
int p[10005],c[10005];
void divide_prime(int n){
    for(int i=2;i*i<=n;i++){
		if(n%i==0){
	    	p[++m]=i;
	    	while(n%i==0)n/=i,c[m]++;
		}
    }
    if(n>1)p[++m]=n,c[m]=1;
}
int quickpow(int a,long long b){
    int cnt=1;
    while(b){
		if(b&1)cnt=((long long)cnt*a)%mod;
		a=((long long)a*a)%mod;
		b>>=1;
    }
    return cnt%mod;
}
int main(){
    a=read(),b=read();
    divide_prime(a);
    for(int i=1;i<=m;i++){
		if((p[i]-1)%mod==0){
	    	ans=(ans*((long long)b*c[i]+1)%mod)%mod;
	    	continue;
		}
		int x=quickpow(p[i],(long long)b*c[i]+1);
		x=(x-1+mod)%mod;
		int y=quickpow(p[i]-1,mod-2);
		ans=((long long)ans*x*y)%mod;
    }
    printf("%d
",ans%mod);
    return 0;
}