（POJ - 1845）Sumdiv（分治法/费马小定理+快速幂）

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

题意：计算a^b的因子的和
分析：将a进行质因子分解，变为p1,p2,p3…pk,显然题目的答案就是(p1^0+p1^1+p1^2+…+p1^a1*b)(p2^0+p2^1+p2^2+p2^3+…+p2^a2*b)…(pk^0+pk^1+pk^2+pk^3+…+pk^ak*b)。

方法1：利用等比公式求和，但是有的数据比较奇葩，这个题和之前写的逆元的题有所不同 59407 2这组数据，59407%9901=1 我写的费马小定理求逆元 之前的题都能AC，但是写这题的时候WA了n次，可能和9901有关？
然后根据别人写的欧几里得的方法写出来，有地方觉得莫名其妙。。 然后吐槽一下，这数据是不是有毒。。
1.费马小+快速幂

#include<cstdio>
#include<iostream>
#include<cstring>
using namespace std;
typedef long long LL;
const int mod=9901;
const int N=1e4+5;
int p[N];
int pnum[N];
int pcnt;
void is_factor(int A)
{
    memset(pnum,0,sizeof(pnum));
    pcnt=0;
    for(int i=2; i*i<=A; i++)
    {
        if(A%i==0)
        {
            p[++pcnt]=i;
            while(A%i==0)
            {
                pnum[pcnt]++;
                A=A/i;
            }
        }
    }
    if(A>1)
    {
        p[++pcnt]=A;
        pnum[pcnt]=1;
    }
}
int pow_mod(int x,int y)
{
    x=x%mod;
    if(x==0) return mod;这个不能删掉。。。   莫名奇妙，无法解释。。  
    int ans=1;
    while(y)
    {
        if(y&1) ans=ans*x%mod;
        x=x*x%mod;
        y=y>>1;
    }
    return ans;
}
LL fermat(LL a,LL p)///费马小定理求逆元
{
    return pow_mod(a,p-2);
}
void solve(int A,int B)
{
    long long int ans=1;
    int x,y;
    for(int i=1; i<=pcnt; i++)
    {
        if(p[i]%mod==1)///特判。。。
        {
            ans=ans*(pnum[i]*B+1)%mod;
            continue;
        }
        ans=ans*(pow_mod(p[i],pnum[i]*B+1)-1)%mod;
        x=fermat(p[i]-1,mod);
        ans=ans*x%mod;
    }
    cout<<ans<<endl;
}
int main()
{
    int A,B;
    while(scanf("%d%d",&A,&B)!=EOF)
    {
        std::ios::sync_with_stdio(false);
        if(A==0)
        {
            cout<<0<<endl;
            continue;
        }
        if(A==1 || B==0)
        {
            cout<<1<<endl;
            continue;
        }
        is_factor(A);
        solve(A,B);
    }
}

2.分治法
运用分治的思想，pi^0+pi^1+pi^2+pi^3+…+pi^ai*b 可以从中部分成两段看成
pi^0+pi^1+pi^2+pi^3+…+pi^(ai*b/2)
pi^(ai*b/2+1)+pi^(ai*b/2+2)+pi^(ai*b/2+3)+…+pi^(ai*b/2)

合并后变成(1+pi^(ai*b/2+1))*？。
这里的问号要分奇偶讨论，详见代码。
递归求解

#include<cstdio>
#include<cstring>
using namespace std;
typedef long long LL;
const int mod=9901;
LL pow_mod(LL a,LL b)
{
    LL ans=1;
    a%=mod;
    while(b>0)
    {
        if(b&1) ans=ans*a%mod;
        b>>=1;
        a=a*a%mod;
    }
    return ans;
}
LL dfs(LL a,LL n)
{
    if(n==0) return 1;
    if(n&1) return ((1+pow_mod(a,n/2+1))*dfs(a,n/2))%mod;
    else return ((1+pow_mod(a,n/2+1))*dfs(a,n/2-1)%mod+pow_mod(a,n/2))%mod;
}
int main()
{
    int A,B;
    while(~scanf("%d%d",&A,&B))
    {
        LL ans=1;
        for(int i=2; i*i<=A; i++)
        {
            int num=0;
            while(A%i==0)
            {
                A/=i;
                num++;
            }
            ans=ans*dfs(i,num*B)%mod;
        }
        if(A>1) ans=ans*dfs(A,B)%mod;
        printf("%lld\n",ans);
    }
    return 0;
}