POJ 1845 唯一分解定理

#include <iostream>
using namespace std;
typedef unsigned long long llong;
const int mod=9901;
int p[4000],len=0,top[32][2],tlen;
bool s[7072]={true,true};
void split(llong n)
{
    tlen=0;
    int i;
    for(i=0;i<len;i++)
        {
            if(n%p[i]==0)
                {
                    top[tlen][0]=0;
                    top[tlen][1]=p[i];
                    while(n%p[i]==0)
                        {
                            top[tlen][0]++;
                            n/=p[i];
                        }
                    tlen++;
                    if(n==1||(n<=7071&&!s[n]))break;
                }
        }
    if(n!=1){top[tlen][0]=1;top[tlen++][1]=n;}
}
llong pow_mod(llong a,llong b,llong c)
{
    if(!b) return 1;
    return pow_mod(a*a%c,b/2,c)*((b&1)?a:1)%c;
}
llong sum(llong a,llong n)
{
   if(n==1) return 1;
   llong k=sum(a,n/2)%mod;
   if(n%2==0) return (k+(pow_mod(a,n/2,mod)*k)%mod)%mod;
   else return (pow_mod(a,n/2,mod)+(k+(pow_mod(a,n/2+1,mod)*k)%mod)%mod)%mod;
}
int main()
{
    llong a,b,ans;
    int i,j;
    for(i=2;i*i<=7071;i++)if(!s[i])for(j=i*i;j<=7071;j+=i)if(!s[j])s[j]=true;
    for(i=2;i<=7071;i++)if(!s[i])p[len++]=i;
    while(cin>>a>>b)
        {
            ans=1;
            split(a);
            for(i=0;i<tlen;i++)
                {
                    ans*=sum(top[i][1],top[i][0]*b+1);
                    ans%=9901;
                }
            cout<<ans<<endl;
        }
    return 0;
}