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

题意：

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

分析：

首先，要先明确一个定理。

整数唯一分解定理：任意大于等于2的正整数都有且只有一种方式写出其质因子的乘积表达式。

A ＝pa11p1a1pa22p2a2pa33p3a3…pannpnan(pi是素数且pi<=pj)

然后A的所有约数和就等于

sum=(p01p10+p11p11+p21p12+…+pa11p1a1)+(p02p20+p12p21+p22p22+…+pa12p2a1)+(p03p30+p13p31+p23p32+…+pa13p3a1)+(p04p40+p14p41+p24p42+…+pa14p4a1)

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

式子写出来了,接下来解决怎么求.

上边的式子每个括号的格式都是一样的，所以会求一个就会求全部了。

很明显每个括号都是一个等比数列。

第一想法自然是等比数列的求和公式。

然鹅，运算过程中要取模的，所以不能有除法，所以求和公式不能用。（不知道逆元能不能求，有兴趣可以试一下）

因此要用二分递归求：

假设现在想要求(p0p0+p1p1+p2p2+…+pnpn)

1.当n是奇数时，括号内有n+1项，有偶数项，把第1项和第（（n+1）/2）项组合，把第2项和第（（n+1）/2+1）项组合，把第3项和第（（n+1）/2+2）项组合…..

可以得到(1+pn/2+11+pn/2+1)××p0p0+(1+pn/2+11+pn/2+1)××p1p1+(1+pn/2+11+pn/2+1)××p2p2+….+(1+pn/2+11+pn/2+1)××pn/2pn/2

再提取公因式(1+pn/2+11+pn/2+1)××(p0p0+p1p1+p2p2+….+pn/2pn/2)

括号里边又是一个等比数列 递归吧

2.当n是偶数时，括号内有n+1项，有奇数项，把第1项和第（n/2+1）项组合，把第2项和第（n/2+2）项组合，把第3项和第（n/2+3）项组合…..

可以得到(1+pn/2+11+pn/2+1)××p0p0+(1+pn/2+11+pn/2+1)××p1p1+(1+pn/2+11+pn/2+1)××p2p2+….+(1+pn/2+11+pn/2+1)××pn/2pn/2

再提取公因式(1+pn/2+11+pn/2+1) ×× (p0p0+p1p1+p2p2+….+pn/2−1pn/2−1)

括号里边又是一个等比数列

emmm…. 简单点说就是 1.当n是奇数时,把式子从中间一切两半，前一半的第1项和后一半的第1项组合，前一半的第2项和后一半的第2项组合…前一半的第n项和后一半的第n项组合。

1.当n是偶数时,先把中间那一项拿掉，然后把式子从中间一切两半，前一半的第1项和后一半的第1项组合，前一半的第2项和后一半的第2项组合…前一半的第n项和后一半的第n项组合。

组合后再提取公因式就行了。

递归终点是n==0的时候返回1。

ok，接下来就简单了。一个是求和的时候要用到快速幂，很基础的一个算法。

快速幂(会的跳过)

然后就是怎么把一个数分解式写出来.

分解质因数，短除法(会的跳过)

分解质因数要用到筛法.

筛法求素数&&线性筛法求素数

因为要求得是NMNM的所有约数，把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;
}