POJ 1845 筛法+分解质因数+快速幂+二分递归求等比数列和

题意：
给你两个数N,M；求$N^M$的所有约数和对9901取模后的结果。(0<=N,M,<=50000000)

分析：

首先，要先明确一个定理。
整数唯一分解定理：任意大于等于2的正整数都有且只有一种方式写出其质因子的乘积表达式。
A ＝$p_1 ^ {a_1}$$p_2 ^ {a_2}$$p_3 ^ {a_3}$…$p_n ^ {a_n}$(pi是素数且pi<=pj)

然后A的所有约数和就等于
sum=($p_1 ^ 0$+$p_1 ^ 1$+$p_1 ^ 2$+…+$p_1 ^ {a_1}$) * ($p_2 ^ 0$+$p_2 ^ 1$+$p_2 ^ 2$+…+$p_2 ^ {a_2}$) * ($p_3 ^ 0$+$p_3^ 1$+$p_3 ^ 2$+…+$p_3 ^ {a_3}$) * ($p_4 ^ 0$+$p_4 ^ 1$+$p_4 ^ 2$+…+$p_4 ^ {a_4}$)…

这个式子的证明：约数和公式证明

式子写出来了,接下来解决怎么求.
上边的式子每个括号的格式都是一样的，所以会求一个就会求全部了。
很明显每个括号都是一个等比数列。
第一想法自然是等比数列的求和公式。
然鹅，运算过程中要取模的，所以不能有除法，所以求和公式不能用。（不知道逆元能不能求，有兴趣可以试一下）
因此要用二分递归求：
假设现在想要求($p ^ 0$+$p ^ 1$+$p ^ 2$+…+$p ^ n$)
1.当n是奇数时，括号内有n+1项，有偶数项，把第1项和第（（n+1）/2）项组合，把第2项和第（（n+1）/2+1）项组合，把第3项和第（（n+1）/2+2）项组合…..
可以得到($1+p ^ {n/2+1}$)$\times$$p ^ 0$+($1+p ^ {n/2+1}$)$\times$$p ^ 1$+($1+p ^ {n/2+1}$)$\times$$p ^ 2$+….+($1+p ^ {n/2+1}$)$\times$$p ^ {n/2}$
再提取公因式($1+p ^ {n/2+1}$)$\times$($p ^ 0$+$p ^ 1$+$p ^ 2$+….+$p ^ {n/2}$)
括号里边又是一个等比数列 递归吧
2.当n是偶数时，括号内有n+1项，有奇数项，把第1项和第（n/2+1）项组合，把第2项和第（n/2+2）项组合，把第3项和第（n/2+3）项组合…..
可以得到($1+p ^ {n/2+1}$)$\times$$p ^ 0$+($1+p ^ {n/2+1}$)$\times$$p ^ 1$+($1+p ^ {n/2+1}$)$\times$$p ^ 2$+….+($1+p ^ {n/2+1}$)$\times$$p ^ {n/2}$
再提取公因式($1+p ^ {n/2+1}$) $\times$ ($p ^ 0$+$p ^ 1$+$p ^ 2$+….+$p ^ {n/2-1}$)
括号里边又是一个等比数列
emmm…. 简单点说就是 1.当n是奇数时,把式子从中间一切两半，前一半的第1项和后一半的第1项组合，前一半的第2项和后一半的第2项组合…前一半的第n项和后一半的第n项组合。
1.当n是偶数时,先把中间那一项拿掉，然后把式子从中间一切两半，前一半的第1项和后一半的第1项组合，前一半的第2项和后一半的第2项组合…前一半的第n项和后一半的第n项组合。
组合后再提取公因式就行了。
递归终点是n==0的时候返回1。
ok，接下来就简单了。一个是求和的时候要用到快速幂，很基础的一个算法。

快速幂(会的跳过)

然后就是怎么把一个数分解式写出来.
分解质因数，短除法(会的跳过)

分解质因数要用到筛法.
筛法求素数&&线性筛法求素数
因为要求得是$N^M$的所有约数，把N的所有约数求出来，再把每个约数的个数乘以M就可以了

还有一点就是这个题的范围是50000000，但是筛法只筛到10000就可以了。
因为你把n所有小于10000的因子全部除掉以后，剩下的n要么是1，要么是一个素数。
因为n最大是50000000，它不可能有两个以上的大于10000的因子。

最后一点，几乎所有运算过程都伴随着取模。
同余定理：
（a+b)%c = (a%c+b%c)%c
(a*b)%c = ((a%c)*(b%c))%c

代码：

#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<iostream>
using namespace std;
#define ll long long
#define Max 10010
#define Mod 9901
bool ispri[Max];
int pri[Max];
int fat[Max][2];
int fatCnt;
ll N,M;

ll _pow(ll a,ll b){//快速幂
    ll res=1;
    while(b){
        if(b&1){
            res=(res*a)%Mod;
        }
        a=(a*a)%Mod;
        b>>=1;
    }
    return res;
}

void isprime(){//线性筛法
    memset(ispri,0,sizeof(ispri));
    memset(pri,0,sizeof(pri));
    ispri[1]=1;
    for(int i=2;i<=Max;i++){
        if(!ispri[i]){
            pri[++pri[0]]=i;
        }
        for(int j=1;j<=pri[0]&&i*pri[j]<=Max;j++){
            ispri[i*pri[j]]=1;
            if(i%pri[j]==0)
                break;
        }
    }
    return;
}
void fen(ll n){//分解质因数
    memset(fat,0,sizeof(fat));
    fatCnt=1;
    for(int i=1;pri[i]*pri[i]<=n;i++){
        while(n%pri[i]==0){
            fat[fatCnt][0]=pri[i];
            fat[fatCnt][1]++;
            n/=pri[i];
        }
        fatCnt++;
    }
    if(n!=1){
        fat[fatCnt][0]=n;
        fat[fatCnt][1]++;
        fatCnt++;
    }
    for(int i=1;i<fatCnt;i++){
        fat[i][1]*=M;
    }
    return;
}                                                                                                                                              
ll sum(ll p,ll n){//二分递归求等比数列和
    if(n==0)return 1;
    if(n&1)return  ((1+_pow(p,n/2+1)%Mod)%Mod*sum(p,n/2)%Mod)%Mod;
    else   return ((1+_pow(p,n/2+1)%Mod)%Mod*sum(p,n/2-1)%Mod+_pow(p,n/2)%Mod)%Mod;
}                    

int main()
{
    isprime();
    while(cin>>N>>M){
        fen(N);
        ll ans=1;
        for(int i=1;i<fatCnt;i++){
            //cout<<fat[i][0]<<" "<<fat[i][1]<<endl;
            ans*=sum(fat[i][0],fat[i][1])%Mod;
            ans%=Mod;
        }
        cout<<ans%Mod<<endl;
    }
    return 0;
}