Sumdiv（poj1845）

题意：求A^B的因子的和。

/*
  首先将A分解 A=p1^a1*p2^a2*...*pn*an
  A^B=p1^a1B*p2^a2B*...*pn*anB
  因子之和sum=(1+p1+p1^2+...+p1^a1B)*...*(1+pn+pn^2+...+pn*anB)
  套用等比数列的公式，再用逆元搞一下。
  求逆元：(a/b)mod m=amod(bm)/b
*/
#include<cstdio>
#include<iostream>
#define N 10010
#define mod 9901
#define lon long long
using namespace std;
int prime[N],f[N],num;
void get_prime(){
    for(int i=2;i<N;i++){
        if(!f[i]) prime[++num]=i;
        for(int j=1;j<=num;j++){
            if(i*prime[j]>=N) break;
            f[i*prime[j]]=1;
            if(i%prime[j]) break;
        }
    }
}
lon msm(lon A,lon B,lon MOD){
    lon base=A,r=0;
    while(B){
        if(B&1) r=(r+base)%MOD;
        base=(base+base)%MOD;
        B>>=1;
    }
    return r;
}
lon ksm(lon A,lon B,lon MOD){
    lon base=A,r=1;
    while(B){
        if(B&1) r=msm(r,base,MOD);
        base=msm(base,base,MOD);
        B>>=1;
    }
    return r;
}
void solve(lon A,lon B){
    lon ans=1;
    for(int i=1;prime[i]*prime[i]<=A;i++){
        if(A%prime[i]==0){
            int num=0;
            while(A%prime[i]==0){
                num++;
                A/=prime[i];
            }
            lon M=(prime[i]-1)*mod;
            ans*=(ksm(prime[i],num*B+1,M)+M-1)/(prime[i]-1);
            ans%=mod;
        }
    }
    if(A>1){
        lon M=mod*(A-1);
        ans*=(ksm(A,B+1,M)+M-1)/(A-1);
        ans%=mod;
    }
    cout<<ans<<endl;
}
int main(){
    lon A,B;
    get_prime();
    while(cin>>A>>B)
        solve(A,B);
    return 0;
}