POJ1845 Sumdiv A^B的约数和（逆元解法或二分乘法）

题目链接：http://poj.org/problem?id=1845

题目大意：计算A^B的所有约数的和对9901取模后的结果。

分析：我们知道，对于一个正整数n，我们有n=(p1)^a1*(p2)^a2*...*(pk)^ak，定义约数和函数σ(n)=∏（pi^(ai+1)-1）/（pi-1）;那么对于A^B，我们有A^B=(p1)^(B*a1)*(p2)^(B*a2)*...*(pk)^(B*ak)，其约数和σ(A^B)=∏（pi^(Ai+B+1)-1）/（pi-1）。对于分母的取模，我们有ans≡a/b(mod m) ==> ans≡a(mod mb)/b。

实现代码如下：

#include <cstdio>
#include <cstring>
using namespace std;
typedef long long LL;
const int MOD=9901;
const int N=10005;
int p[N],cnt;
bool isp[N];
void Init()
{
    cnt=0;
    memset(isp,true,sizeof(isp));
    for(int i=2;i<N;i++)
      if(isp[i])
      {
          p[cnt++]=i;
          for(int j=i+i;j<N;j+=i)
            isp[j]=false;
      }
}
LL multi(LL a,LL b,LL m)
{
    LL ans=0;
    a%=m;
    while(b)
    {
        if(b&1) ans=(ans+a)%m;
        b>>=1;
        a=(a+a)%m;
    }
    return ans;
}
LL quick_mod(LL a,LL b,LL m)
{
    LL ans=1;
    a%=m;
    while(b)
    {
        if(b&1) ans=multi(ans,a,m);
        b>>=1;
        a=multi(a,a,m);
    }
    return ans;
}
void Solve(LL A,LL B)
{
    LL ans=1;
    for(int i=0;i<cnt&&p[i]*p[i]<=A;i++)
    {
        if(A%p[i]==0)
        {
            int num=0;
            while(A%p[i]==0)
            {
                num++;
                A/=p[i];
            }
            LL M=(p[i]-1)*MOD;
            ans*=(quick_mod(p[i],num*B+1,M)-1)/(p[i]-1);
            ans%=MOD;
        }
    }
    if(A>1)
    {
        LL M=(A-1)*MOD;
        ans*=(quick_mod(A,B+1,M)-1)/(A-1);
        ans%=MOD;
    }
    printf("%I64d\n",ans);
}
int main()
{
    LL A,B;
    Init();
    while(scanf("%I64d%I64d",&A,&B)!=-1)
      Solve(A,B);
    return 0;
}

对于等比数列的求和，我们还可以用二分乘法来做，具体做法可以参考：点击打开链接。

假设有等比数列1,a,a^2,...,a^n:

（1）如果n为奇数，那么(n+1)为偶数，此时我们有:

S(n)=1+a+a^2+...+a^n=1+a+a^2+...+a^((n-1)/2)+a^((n-1)/2+1)+...+a^((n-1)/2+(n-1)/2)+a^((n-1)/2+(n-1)/2+1)=(1+a^((n-1)/2+1) * ( 1+a+a^2+...+a^((n-1)/2) )=(1+a^((n-1)/2+1)*S((n-1)/2);

（2）如果n为偶数，那么(n+1)为奇数，此时我们有：

S(n)=1+a+a^2+...+a^n=1+a+a^2+...+a^(n/2-1)+a^(n/2)+a^(n/2+1)+...+a^(n/2+n/2)=(1+a^(n/2+1) * ( 1+a+a^2+...+a^(n/2-1) )=(1+a^(n/2+1)*S(n/2-1);

实现代码如下：

#include <cstdio>
#include <cstring>
using namespace std;
typedef long long LL;
const int MOD=9901;
const int N=10005;
int p[N],cnt;
bool isp[N];
void Init()
{
    cnt=0;
    memset(isp,true,sizeof(isp));
    for(int i=2;i<N;i++)
      if(isp[i])
      {
          p[cnt++]=i;
          for(int j=i+i;j<N;j+=i)
            isp[j]=false;
      }
}
LL quick_mod(LL a,LL b,LL m)
{
    LL ans=1;
    a%=m;
    while(b)
    {
        if(b&1) ans=ans*a%m;
        b>>=1;
        a=a*a%m;
    }
    return ans;
}
LL Sum(LL a,LL n)
{
    if(n==0) return 1;
    if(n&1)
    {
        return ( (1+quick_mod(a,(n-1)/2+1,MOD) ) * Sum(a,(n-1)/2)%MOD)%MOD;
    }
    else
    {
        return ( (1+quick_mod(a,n/2+1,MOD)) * Sum(a,n/2-1)%MOD + quick_mod(a,n/2,MOD) )%MOD;
    }
}
void Solve(LL A,LL B)
{
    LL ans=1;
    for(int i=0;i<cnt&&p[i]*p[i]<=A;i++)
    {
        if(A%p[i]==0)
        {
            int num=0;
            while(A%p[i]==0)
            {
                num++;
                A/=p[i];
            }
            ans*=Sum(p[i],num*B)%MOD;
            ans%=MOD;
        }
    }
    if(A>1)
    {
        ans*=Sum(A,B)%MOD;
        ans%=MOD;
    }
    printf("%I64d\n",ans);
}
int main()
{
    LL A,B;
    Init();
    while(scanf("%I64d%I64d",&A,&B)!=-1)
      Solve(A,B);
    return 0;
}