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


评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值