BZOJ3994: [SDOI2015]约数个数和

BZOJ3994: [SDOI2015]约数个数和


题目描述

传送门

题目分析

求的东西简明扼要,
\[\sum_{i=1}^{N}\sum_{j=1}^{m}d(ij)\]

但是有个需要知道的是
\[d(ij)=\sum_{x\mid i}\sum_{y\mid j}[gcd(x,y)=1]\]
然后可以开始开心的推式子。

根据反演套路,设两个函数
\[ f(d)=\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=d] \]
\[ F(n)=\sum_{n\mid d}f(d)=\lfloor \frac Nd \rfloor \lfloor \frac Md \rfloor \]
根据莫比乌斯反演可以得到:
\[ f(n)=\sum_{n\mid d}\mu(\frac dn)F(d) \]
\(d(i,j)\)代回原来的式子并化简:
\[ \begin{align} Ans&=\sum_{i=1}^{N}\sum_{j=1}^{M}d(ij) \\ &=\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{x\mid i}\sum_{y\mid j}[gcd(x,y)=1]\ 根据\mu的性质把它代进去\\ &=\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{x\mid i}\sum_{y\mid j}\sum_{d\mid gcd(x,y)}\mu(d)\ 然后更换枚举约数为枚举d \\ &=\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{x\mid i}\sum_{y\mid j}\sum_{d=1}^{min(N,M)}\mu(d)\times [d\mid gcd(x,y)]\\ &=\sum_{d=1}^{min(N,M)}\mu(d)\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{x\mid i}\sum_{y\mid j}[d\mid gcd(x,y)]\\ &=\sum_{d=1}^{min(N,M)}\mu(d)\sum_{x=1}^{N}\sum_{y=1}^{M}[d\mid gcd(x,y)]\lfloor \frac {N}{x} \rfloor \lfloor \frac{M}{y}\rfloor \\ &=\sum_{d=1}^{min(N,M)}\mu(d)\sum_{x=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{M}{y}\rfloor}\lfloor\frac{N}{dx}\rfloor\lfloor\frac{M}{dy}\rfloor\\ &=\sum_{d=1}^{min(N,M)}\mu(d)(\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{dx}\rfloor)(\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\lfloor\frac{m}{dy}\rfloor) \end{align} \]
明显,推到这里了之后就可以在\(O(n)\)的时间内求出了。

对于多组数据,可以进行整除分块,再将这个式子优化成\(O(\sqrt n)\)的。

是代码呢

#include <bits/stdc++.h>
using namespace std;
const int MAXN=1e5+7;
#define ll long long
ll sum[MAXN],g[MAXN];
int mu[MAXN],prime[MAXN];
bool vis[MAXN];
int n,m;
inline void get_mu(int N)
{
    mu[1]=1;
    for(int i=2;i<=N;i++){
        if(!vis[i]){
            mu[i]=-1;
            prime[++prime[0]]=i;
        }
        for(int j=1;j<=prime[0];j++){
            if(prime[j]*i>N) break;
            vis[prime[j]*i]=1;
            if(i%prime[j]==0) break;
            else mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=N;i++) sum[i]=sum[i-1]+mu[i];
    for(int i=1;i<=N;i++){
        ll ans=0;
        for(int l=1,r;l<=i;l=r+1){
            r=(i/(i/l));
            ans+=1ll*(r-l+1)*(i/l);
        }
        g[i]=ans;
    }
}
inline int read()
{
    int x=0,c=1;
    char ch=' ';
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    while(ch=='-')c*=-1,ch=getchar();
    while(ch>='0'&&ch<='9')x=x*10+ch-'0',ch=getchar();
    return x*c;
}
int main()
{
    int T=read();
    get_mu(50000);
    while(T--){
        n=read(),m=read();
        int r=0,mx=min(n,m);
        ll ans=0;
        for(int l=1;l<=mx;l=r+1){
            r=min((n/(n/l)),(m/(m/l)));
            ans+=1ll*(sum[r]-sum[l-1])*(g[n/l])*g[m/l];
        }
        printf("%lld\n", ans);
    }
}

转载于:https://www.cnblogs.com/victorique/p/10385823.html

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值