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

题目：传送门

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

分析:
这道题折腾了很久啊，值得写一写报告。
首先一个大于$1$的正整数$X$有唯一分解：

$$X=p_1^{e_1}*p_2^{e_2}*...*p_n^{e_n}，其中p_n为素数。$$
那么 $X$的所有因子之和为：
$$sum=(1+p_1^1+p_1^2+...+p_1^{e_1})*(1+p_2^1+p_2^2+...+p_2^{e_2})*...*(1+p_n^1+p_n^2+...+p_n^{e_n})$$
所以我们先求 $A$的素因子分解，然后答案就是：
$$ans=(1+p_1^1+p_1^2+...+p_1^{e_1*B})*(1+p_2^1+p_2^2+...+p_2^{e_2*B})*...*(1+p_n^1+p_n^2+...+p_n^{e_n*B})\%mod$$
求多项式 $(1+x^1+x^2+...+x^n)\%mod$有两种方法：
1）递归二分：
当 $n$为奇数时， $$sum\%mod=(1+x^1+x^2+...+x^{n/2})*(1+x^{n/2+1})\%mod$$ 可以看到，前面是原式的一半，可以递归求解，而后面是指数式，可以用快速幂算法。
当 $n$为偶数时， $$sum\%mod=((1+x^1+x^2+...+x^{n/2-1})*(1+x^{n/2+1})+x^{n/2})\%mod$$ 同样使用递归和快速幂求解。
2）等比通项：
显然 $$sum\%mod=\frac {x^{n+1}-1}{x-1}\%mod$$ 形如 $\frac{a}{b}\%mod$的分式取模也有两种方法，一种是求分母 $b$在模$mod$意义下的逆元，那么 $\frac{a}{b}\%mod=a*b^{-1}\%mod$，这里因为 $mod=9901$是素数，所以由欧拉定理： $b^{\phi(mod)}\equiv1\%mod$可知 $b^{-1}\equiv b^{mod-2}\%mod$另一种是变换模值，即 $\frac{a}{b}\%mod=\frac{a\%(mod*b)}{b}\%mod$，在本题中，有可能某些素因子 $p$满足$(p-1)\%mod==0$所以此时 $p-1$不存在模 $mod$意义下的逆元，此时只能使用变换模值法，而使用变换模值法时要注意 $mod*b$有可能很大，在后续计算时会造成溢出错误，这里我们采用 $\%mod$逐位连乘的技巧，具体看代码吧。

ps:在一篇题解中，大神给出了求$1 \to p(p为奇素数)$模$p$的所有逆元的一个递推式，用这个递推式就可以在$O(p)$复杂度里求逆元了。这对于 p <script type="math/tex" id="MathJax-Element-165">p</script>很大时是很有用的。

代码：

/*等比通项*/
#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;
}