POJ 1845 Sumdiv <数论（逆元 / 二分递归）>

题目：传送门

题目大意：输入两个数A，B求AB的所有因子之和。

分析:

这道题折腾了很久啊，值得写一写报告。

首先一个大于1的正整数X有唯一分解：

X=pe11∗pe22∗...∗penn，其中pn为素数。

那么X的所有因子之和为：

sum=(1+p11+p21+...+pe11)∗(1+p12+p22+...+pe22)∗...∗(1+p1n+p2n+...+penn)

所以我们先求A的素因子分解，然后答案就是：

ans=(1+p11+p21+...+pe1∗B1)∗(1+p12+p22+...+pe2∗B2)∗...∗(1+p1n+p2n+...+pen∗Bn)%mod

求多项式(1+x1+x2+...+xn)%mod有两种方法：

1）递归二分：

当n为奇数时，sum%mod=(1+x1+x2+...+xn/2)∗(1+xn/2+1)%mod可以看到，前面是原式的一半，可以递归求解，而后面是指数式，可以用快速幂算法。

当n为偶数时，sum%mod=((1+x1+x2+...+xn/2−1)∗(1+xn/2+1)+xn/2)%mod同样使用递归和快速幂求解。

2）等比通项：

显然sum%mod=xn+1−1x−1%mod形如ab%mod的分式取模也有两种方法，一种是求分母b在模mod意义下的逆元，那么ab%mod=a∗b−1%mod，这里因为mod=9901是素数，所以由欧拉定理：bϕ(mod)≡1%mod可知b−1≡bmod−2%mod另一种是变换模值，即ab%mod=a%(mod∗b)b%mod，在本题中，有可能某些素因子p满足(p−1)%mod==0所以此时p−1不存在模mod意义下的逆元，此时只能使用变换模值法，而使用变换模值法时要注意mod∗b有可能很大，在后续计算时会造成溢出错误，这里我们采用%mod逐位连乘的技巧，具体看代码吧。

ps:在一篇题解中，大神给出了求1→p(p为奇素数)模p的所有逆元的一个递推式，用这个递推式就可以在O(p)复杂度里求逆元了。这对于p很大时是很有用的。

代码：

/*等比通项*/
#include <iostream>
#include <cstring>
#include <string>
using namespace std;

typedef long long LL;
LL a,b;
LL p,e;
LL m=9901;

LL multi(LL a,LL b,LL mod){
LL res=0;
a%=mod;
b%=mod;
while(b){
if(b&1) res=(res+a)%mod;
a=(a*2)%mod;
b>>=1;
}
return res;
}

LL power(LL a,LL b,LL mod){
LL res=1;
a%=mod;
while(b){
if(b&1) res=multi(res,a,mod);
a=multi(a,a,mod);
b>>=1;
}
return res;
}

LL sum(LL p,LL n, mod){
if(!n) return 1;
else{
if(n&1) return sum(p,n/2)*(power(p,n/2+1,mod)+1)%mod;
else return sum(p,n/2-1)*(power(p,n/2+1,mod)+1)*power(p,n/2,mod)%mod;
}
}

void solve(){
LL ans=1;
for(int i=2;i*i<=a;){
if(a%i==0){
p=i;
e=0;
while(!(a%i)){
e++;
a/=i;
}
if((p-1)%m==0){
LL M=m*(p-1);
ans=ans*(power(p,e*b+1,M)+M-1)%M/(p-1)%m;
}
else ans=ans*(power(p,e*b+1,m)-1)*power(p-1,m-2,m)%m;
}
if(i==2) i++;
else i+=2;
}
if(a>1){
p=a; e=1;
if((p-1)%m==0){
LL M=m*(p-1);
ans*=(power(p,e*b+1,M)+M-1)%M/(p-1)%m;
}
else ans=ans*(power(p,e*b+1,m)-1)*power(p-1,m-2,m)%m;
}
cout<<(ans+m)%m<<endl;
}

int main(){
while(cin>>a>>b){
solve();
}
return 0;
}

/*递归二分*/
#include <iostream>
#include <cstring>
#include <string>
using namespace std;

typedef long long LL;
LL a,b;
LL p,e;
LL mod=9901;

LL power(LL a,LL b){
LL res=1;
a%=mod;
while(b){
if(b&1) res=res*a%mod;
a=a*a%mod;
b>>=1;
}
return res;
}

LL sum(LL p,LL n){
if(!n) return 1;
else{
if(n&1) return sum(p,n/2)*(power(p,n/2+1)+1)%mod;
else return (sum(p,n/2-1)*(power(p,n/2+1)+1)+power(p,n/2))%mod;
}
}

void solve(){
LL ans=1;
for(int i=2;i*i<=a;){
if(a%i==0){
p=i;
e=0;
while(!(a%i)){
e++;
a/=i;
}
ans=ans*sum(p,e*b)%mod;
}
if(i==2) i++;
else i+=2;
}
if(a>1){
p=a; e=1;
ans=ans*sum(p,e*b)%mod;
}
cout<<(ans+mod)%mod<<endl;
}

int main(){
while(cin>>a>>b){
solve();
}
return 0;
}