poj 1845 Sumdiv

DescriptionConsider two natural numbers A and B. Let S be the sum of all natural divisors of A^B. Determine S modulo 9901 (the rest of the division of S by 9901).InputThe only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.OutputThe only line of the output will contain S modulo 9901.Sample Input2 3Sample Output15Hint2^3 = 8. The natural divisors of 8 are: 1,2,4,8. Their sum is 15. 15 modulo 9901 is 15 (that should be output). 

#include<cstdio>#include<cstdlib>#include<cstring>#include<cmath>#include<iostream>#include<algorithm>#include<vector>#include<map>#include<set>#include<queue>#include<stack>using namespace std;const int inf=0x3f3f3f3f;const double pi=3.14159265358979323846264338327950288L;const double eps=1e-6;#define MOD 9901int prime[100000+100]={0},num_prime=0;void get_prime(int m)  //求出m以内的质数并保存到prime数组中{    long long abc[100000]={1,1},i,j;    for(i=2;i<m;i++)    {        if(abc[i]==0)            prime[num_prime++]=i;        for(j=0;j<num_prime&&i*prime[j]<m;j++)        {            abc[i*prime[j]]=1;            if((i%prime[j])==0)                break;        }    }}long long qpow(long long a,long long n)  //求a^n，对MAX取模{    long long ans=1;    if(n==0)        return 1;    a=a%MOD;    while(n>0)    {        if(n%2==1)            ans=(ans*a)%MOD;        n/=2;        a=(a*a)%MOD;    }    return ans;}long long qsum(long long a,long long n){    if(n==0)        return 1;    long long cur,ti=qsum(a,(n-1)/2);    if(n&1)    {        cur = qpow(a,(n+1)/2);        ti=(ti+ti*cur)%MOD;    }    else    {        cur=qpow(a,(n+1)/2);        ti=(ti+ti*cur+qpow(a,n))%MOD;    }    return ti;}int main(){    int i,j,k;    long long a,b,ans;    get_prime(10000);    while(cin>>a>>b)    {        ans=1;        for(i=0;i<num_prime;i++)        {            if(a%prime[i]==0)            {                j=0;                while(a%prime[i]==0)                {                    j++;                    a/=prime[i];                }                ans=(ans*qsum(prime[i],j*b))%MOD;            }        }        if(a>1)        {            ans=(ans*qsum(a,b))%MOD;        }        cout<<ans<<endl;    }    return 0;}