poj1845--Sumdiv(数论篇3--真滴是数论啊。。。。)

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std ;
#define MOD 9901
#define LL __int64
LL vis[1100000] , pri[1100000] , cnt ;
LL a[1100000] , num[1100000] ;
void get_pri()
{
    LL i , j ;
    memset(vis,0,sizeof(vis)) ;
    cnt = 0 ;
    for(i = 2 ; i < 1010000 ; i++)
    {
        if( !vis[i] )
            pri[cnt++] = i ;
        for(j = 0 ; j < cnt ; j++)
        {
            if( pri[j] * i > 1010000 ) break ;
            vis[ pri[j]*i ] = 1 ;
            if(i%pri[j]==0) break ;
        }
    }
    return ;
}
LL pow(LL a,LL k,LL mod)
{
    if( k == 0 ) return 1 ;
    LL b = pow(a,k/2,mod) ;
    b = b*b % mod ;
    if( k%2 ) b = b*a%mod ;
    return b ;
}
int main()
{
    get_pri() ;
    LL n , m , i , j ;
    LL ans , temp ;
    while( scanf("%I64d %I64d", &n, &m) != EOF )
    {
        ans = 1 ;
        memset(num,0,sizeof(num)) ;
        j = 0 ;
        for(i = 0 ; i < cnt ; i++)
        {
            while( n % pri[i] == 0 )
            {
                n /= pri[i] ;
                a[j] = pri[i] ;
                num[j]++ ;
            }
            if( num[j] )
                j++ ;
            if( n == 1 ) break ;
        }
        if( n > 1 )
        {
            a[j] = n ;
            num[j++] = 1 ;
        }
        for(i = 0 ; i < j ; i++)
        {
            if( (a[i]-1)%MOD == 0 )
            {
                temp = pow(a[i],num[i]*m+1,MOD*(a[i]-1)) - 1 ;
                if( temp < 0 ) temp += MOD*(a[i]-1) ;
                ans = ans * ( temp / (a[i]-1) % MOD ) % MOD ;
            }
            else
            {
                temp = pow(a[i],num[i]*m+1,MOD) - 1 ;
                if( temp < 0 ) temp += MOD ;
                ans = ans * (temp) % MOD * pow(a[i]-1,MOD-2,MOD) % MOD ;

            }
        }
        printf("%I64d\n", ans) ;
    }
    return 0 ;
}