POJ 1845 Sumdiv

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

题意：

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

思路：

记一个数的因子和为sigma(x)，如果x为合数，将x进行质因子分解
x=$p_1^{a_1}*p_2^{a_2}*……*p_n^{a_n}$
sigma(x)=$sigma(p_1^{a_1})*sigma(p_2^{a_2})*……*sigma(p_n^{a_n})$
=($1+p_1+p_1^2+...+p_1^{a_1}) * …… * (1+p_n+p_n^2+...+p_n^{a_n}$)
对于每个因子求$1+p_n+p_n^2+...+p_n^{a_n}$再相乘即可
对于这个式子的求法有两种方法：
①递归二分法，具体看代码版本三的sum函数(个人认为这个方法最好，不涉及逆元什么的，简单)
②利用等比数列的求和公式可以求得$(p^{a_n+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;
}