poj 1845 Sumdiv

#include<iostream>        //求 A^B 的所有约数之和
using namespace std;
const __int64 mod=9901;
__int64 fac[1000],num[1000],rear;
__int64 root(__int64 a,__int64 b)    //二分法求解 a^b 
{
    if(b==0)
        return 1;
    else if(b==1)
        return a%mod;
    __int64 c=root(a,b/2),d=c*c%mod;
    if(b%2==0)
        return d;
    else
        return d*a%mod;
}
void factor(__int64 a)    //分解质因子
{
    __int64 i=0;
    for(__int64 j=2;j*j<=a;++j)
        if(a%j==0)
        {
            fac[i]=j;
            while(a%j==0)
            {
                a=a/j;
                num[i]++;
            }
            i++;
        }
    if(a!=1)
        num[i]=1,fac[i]=a,i++;
    rear=i;
}

//1+p+p^2+...p^n=(1+p+p^2+...+p^(n/2))*(1+p^(n/2+1))   n为奇数
//1+p+p^2+...p^n=(1+p+p^2+...+p^(n/2))*(1+p^(n/2+1))-p^(n+1)  n为偶数
__int64 factorial(__int64 p,__int64 n)        //求解 1+p+p^2+...p^n
{
    if(n==0)
        return 1;
    if(n%2==1)
        return factorial(p,n/2)*(1+root(p,n/2+1))%mod;
    else
        return (factorial(p,n/2)*(1+root(p,n/2+1))-root(p,n+1)+mod)%mod;
}
int main()
{
    __int64 a,b,s=1;
    cin>>a>>b;
    factor(a);
    for(__int64 i=0;i<rear;++i)
    {
        s=s*factorial(fac[i],num[i]*b)%mod;
    }
    cout<<s<<endl;
    return 0;
}
//将A进行质因数分解得到p1, p2, p3....数量分别是n1, n2, n3...
//A^B的所有约数的和是
//(p1^0 + p1^1 + ......+p1^(n1 * B)) * (p2 ^0 + p2^1 + ....p2^(n2 * B))*...