poj 1845 求A^B的约数之和

题意：
求A^B的所有约数之和 Mod 9901。

思路：
大数模运算。两个最基本公式：(A*B)%C = ((A%C)(B%C))%C 和 (A+B)%C = ((A%C)+(B%C))%C 。用__int64的原因为 n = cnt[i] B (cnt[i]为A第i个素因子的个数)可能会超int。

1: 对A进行素因子分解得

 A = p1^a1 * p2^a2 * p3^a3 *...* pn^an.

故 A^B = p1^(a1*B) * p2^(a2*B) … pn^(an*B);

2：A^B的所有约数之和为：

 sum = [1+p1+p1^2+...+p1^(a1*B)] * [1+p2+p2^2+...+p2^(a2*B)] *...* [1+pn+pn^2+...+pn^(an*B)].

如 200 = 2^3 * 5^2 : sum(200) = [1 + 2 + 4 + 8] * [1 + 5 + 25].

3: 求等比数列1+pi+pi^2+pi^3+…+pi^n可以由递归形式的二分求得：（模运算不能用等比数列和公式！）
若n为奇数,一共有偶数项，则：

  1 + p + p^2 + p^3 +...+ p^n
  = (1+p^(n/2+1)) + p * (1+p^(n/2+1)) +...+ p^(n/2) * (1+p^(n/2+1))
  = (1 + p + p^2 +...+ p^(n/2)) * (1 + p^(n/2+1))

如：1 + p + p^2 + p^3 + p^4 + p^5 = (1 + p + p^2) * (1 + p^3)
若n为偶数,一共有奇数项,则:

  1 + p + p^2 + p^3 +...+ p^n
  = (1+p^(n/2+1)) + p * (1+p^(n/2+1)) +...+ p^(n/2-1) * (1+p^(n/2+1)) + p^(n/2)
  = (1 + p + p^2 +...+ p^(n/2-1)) * (1+p^(n/2+1)) + p^(n/2);

如：1 + p + p^2 + p^3 + p^4 = (1 + p) * (1 + p^3) + p^2
上述分析copy自：http://blog.sina.com.cn/s/blog_6635898a0100omcn.html

#include<iostream>
#include<cstdio>
#include<set>
using namespace std;
typedef long long ll;
const int mod=9901;
int qmod(ll a,ll n) //快速幂取模
{
    int ans=1;
    for(;n;n>>=1){
        if(n&1)ans=(ll)(ans*a)%mod;
        a=a*a%mod;
    }
    return ans;
}
ll sum(ll p,ll n) //递归求和
{
    if(n==0)return 1;
    if(n&1)return ((1+qmod(p,n/2+1))*sum(p,n/2))%mod;
    else return ((1+qmod(p,n/2+1))*sum(p,n/2-1)+qmod(p,n/2))%mod;
}
int p[100],cnt[100],tot;
void getfact(int a)  //分解质因数
{
    tot=0;
    for(int i=2;i*i<=a;i++){
        if(a%i==0){
            p[tot]=i;
            while(a%i==0){
                cnt[tot]++;
                a/=i;
            }
            tot++;
        }
    }
    if(a!=1){
        p[tot]=a;
        cnt[tot++]=1;
    }
}
int main()
{
    //freopen("f.txt","r",stdin);
    int a,b;
    scanf("%d%d",&a,&b);
    getfact(a);
    int ans=1;
    for(int i=0;i<tot;i++)
        ans=ans*(sum(p[i],(ll)cnt[i]*b)%mod)%mod;
    printf("%d\n",ans);
    return 0;
}