poj 1845 Sumdiv

Description

Consider 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).
Input

The only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.
Output

The only line of the output will contain S modulo 9901.
Sample Input

2 3
Sample Output

15
Hint

2^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 9901

int 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;
}