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;}