POJ 1845 Sumdiv

转载请注明出处，谢谢http://blog.csdn.net/bigtiao097?viewmode=contents

题意：

给两个数A,B，求A^B所有因子的和%9901，A、B的范围是5e7

思路：

记一个数的因子和为sigma(x)，如果x为合数，将x进行质因子分解
x=pa11∗pa22∗……∗pann
sigma(x)=sigma(pa11)∗sigma(pa22)∗……∗sigma(pann)
=(1+p1+p21+...+pa11)∗……∗(1+pn+p2n+...+pann)
对于每个因子求1+pn+p2n+...+pann再相乘即可
对于这个式子的求法有两种方法：
①递归二分法，具体看代码版本三的sum函数(个人认为这个方法最好，不涉及逆元什么的，简单)
②利用等比数列的求和公式可以求得(pan+1−1)/(p−1)
对于这个式子的求法又有两种，一是利用a/b % c = (a % (b*c))/ b ，这样不用求逆元，但是在快速幂的时候可能会爆long long，需要在快速幂中加入快速乘，详见代码版本一；二是求p-1对9901的逆元，但是这里有几种特殊情况，比如p-1是9901的倍数就不能直接去求逆元，需要特殊处理，详见代码版本二

---

具体代码如下：
Result：Accepted    Time : 16ms
版本一

#include<iostream>using namespace std;const int mod = 9901;const int maxn = 5e7+5;typedef long long ll;ll fac[20];//第i个因数ll num[20]={0};//第i个因数的几次幂int cnt=0;//注意不要忘记初始化void divide(int n){    for(int i=2;i*i<=n;)    {         if(n%i == 0)        {            fac[++cnt]=i;            num[cnt]= 0;            while(n%i == 0)            {                num[cnt]++;                n /= i;            }        }        if(i==2)            i++;        else            i+=2;    }    if(n!=1)    {        fac[++cnt] = n;        num[cnt] = 1;    }}ll mod_mul(ll x,ll n,ll p){    x%=p;    ll res = 0;    while(n>0)    {        if(n&1) res = (res+x)%p;        x = (x<<1)%p;        n>>=1;    }    return res;}ll mod_pow(ll x,ll n,ll p){    x%=p;    ll res = 1;    while(n>0)    {        if(n&1) res = mod_mul(res,x,p)%p;        x =mod_mul(x,x,p)%p;        n>>=1;    }    return res;}int main(){    ll a,b;    cin>>a>>b;    if(a==0)    {        cout<<0<<endl;        return 0;    }    if(a==1||b==0)    {        cout<<1<<endl;        return 0;    }    divide(a);    ll ans = 1;    for(int i=1;i<=cnt;i++)    {        ll p = (fac[i]-1)*mod;        ans *= ((mod_pow(fac[i],num[i]*b+1,p)-1+p)%p)/(fac[i]-1);        ans %= mod;    }    cout<<ans<<endl;}

版本二：

#include<iostream>using namespace std;const int mod = 9901;const int maxn = 5e7+5;typedef long long ll;int fac[20];//第i个因数ll num[20]={0};//第i个因数的几次幂int cnt=0;//注意不要忘记初始化void divide(int n){    for(int i=2;i*i<=n;)    {         if(n%i == 0)        {            fac[++cnt]=i;            num[cnt]= 0;            while(n%i == 0)            {                num[cnt]++;                n /= i;            }        }        if(i==2) i++;        else i+=2;    }    if(n!=1)    {        fac[++cnt] = n;        num[cnt] = 1;    }}ll mod_mul(ll x,ll n,ll p){    x%=p;    ll res = 0;    while(n>0)    {        if(n&1) res = (res+x)%p;        x = (x<<1)%p;        n>>=1;    }    return res;}ll mod_pow(ll x,ll n,ll p){    x%=p;    ll res = 1;    while(n>0)    {        if(n&1) res = mod_mul(res,x,p)%p;        x =mod_mul(x,x,p)%p;        n>>=1;    }    return res;}int main(){    ll a,b;    cin>>a>>b;    if(a==0)    {        cout<<0<<endl;        return 0;    }    if(a==1||b==0)    {        cout<<1<<endl;        return 0;    }    divide(a);    ll ans = 1;    for(int i=1;i<=cnt;i++)    {        fac[i]%= mod;        if(fac[i]==1)        {            ans*= num[i]*b+1;            ans%=mod;            continue;        }        if(fac[i]==0)            continue;        ans *= ((mod_pow(fac[i],num[i]*b+1,mod)-1+mod)%mod)*(mod_pow(fac[i]-1,mod-2,mod));        ans %= mod;    }    cout<<ans<<endl;}

版本三：

#include<iostream>using namespace std;const int mod = 9901;typedef long long ll;int fac[20];//第i个因数int num[20];//第i个因数的几次幂int cnt=0;//注意不要忘记初始化void divide(int n){    for(int i=2;i*i<=n;)    {         if(n%i == 0)        {            fac[++cnt]=i;            num[cnt]= 0;            while(n%i == 0)            {                num[cnt]++;                n /= i;            }        }        if(i==2)            i++;        else            i+=2;    }    if(n!=1)    {        fac[++cnt] = n;        num[cnt] = 1;    }}ll mod_pow(ll x,ll n){    x%=mod;    ll res = 1;    while(n>0)    {        if(n&1) res = res*x%mod;        x = x*x%mod;        n>>=1;    }    return res;}ll inv(ll a,ll mod){    return mod_pow(a,mod-2)%mod;}ll sum(ll p,ll n)//计算1+p+p^2+````+p^n{    if(p==0)return 0;    if(n==0)return 1;    if(n&1)//奇数        return ((1+mod_pow(p,n/2+1))%mod*sum(p,n/2)%mod)%mod;    else    return((1+mod_pow(p,n/2+1))%mod*sum(p,n/21)+mod_pow(p,n/2)%mod)%mod;}int main(){    int a,b;    cin>>a>>b;    if(a==0)    {        cout<<0<<endl;        return 0;    }    if(a==1||b==0)    {        cout<<1<<endl;        return 0;    }    ll ans = 1;    divide(a);    for(int i=1;i<=cnt;i++)    {        num[i]*=b;        ans = (ans * sum(fac[i],num[i]))%mod;    }    cout<<ans<<endl;}