POJ 1845 Sumdiv

解题思路:

对A质因子分解A=p1^n1*p2^n2*...*pm^nm,则A^B=p1^(n1*B)*p2^(n2*B)*...*pm^(nm*B);

A^B的因子之和为(1+p1^1+...+p1^n1*B)*(1+p2^1+...+p2^n2*B)*...*(1+pm^1+pm^2+...+pm^nm*B);

1+pi^1+...+pi^m=(1+pi^m-1)/(pi-1),因此因子之和可以化成p/q的形式，且q,M互质时:(p/q)mod(M)=p*q^(-1)mod(M)

q^(-1)为q的乘法逆元，M为质数时q^(-1)=q^(M-2)，所以(p/q)mod(M) = (p * q^(p-2))mod(M)

注:对(pi-1)%M=0的情况需特别处理Line 31~35

#include <iostream>
using namespace std;
#define MAXN 7100
#define MOD 9901
#define MI __int64
MI Cmod(MI a, MI b, MI m)
{
    if(b==0)return 1;
    if(b&1)return a * Cmod(a, b-1, m) % m;
    else return Cmod((a*a)%m, b>>1, m);
}
int main()
{
    bool IsPri[MAXN];
    int a,b,i,j,p,q,pri[1000],m[100],n[100];
    MI num,den,mod,temp,ans;
    memset(IsPri, 0x01, sizeof(IsPri));
    for(IsPri[0]=IsPri[1]=0,i=2;i<85;i++)
        if(IsPri[i])
            for(j=i;i*j<MAXN;j++)IsPri[i*j]=0;
    for(p=0,i=2;i<MAXN;i++)if(IsPri[i])pri[p++]=i;
    scanf("%d %d", &a, &b);
    memset(n,0,sizeof(n));
    for(q=i=0;i<p&&a!=1;i++)
        if(!(a%pri[i]))
            for(m[q++]=pri[i];!(a%pri[i]);a=a/pri[i],n[q-1]++);
    if(a!=1){m[q++]=a,n[q-1]=1;}
    for(num=den=1,i=0;i<q;i++)
    {
        mod = MOD;
        if((m[i]-1)%MOD==0) mod *= (m[i]-1);
        temp=Cmod(m[i],n[i]*b+1,mod);
        num=(num*(temp+mod-1))%mod;
        den = den*(m[i]-1)%mod;
        if((m[i]-1)%MOD==0)num/=MOD,den/=MOD;
    }
    den = Cmod(den, MOD-2,MOD);//逆元的求取
    ans = (den * num)%MOD;
    printf("%I64d\n", ans);
    return 0;
}