#莫比乌斯反演,整除分块,线性筛#bzoj 3994 洛谷 3327 [SDOI2015]约数个数和

题目

d ( x ) d(x) d(x) x x x的约数个数,给定 n , m n,m n,m,求 ∑ i = 1 n ∑ j = 1 m d ( i j ) \sum^n_{i=1}\sum^m_{j=1}d(ij) i=1nj=1md(ij)


分析

又到了推式子的过程了
a n s = ∑ i = 1 n ∑ j = 1 m ∑ x ∣ i ∑ y ∣ j [ g c d ( x , y ) = 1 ] ans=\sum^n_{i=1}\sum^m_{j=1}\sum_{x|i}\sum_{y|j}[gcd(x,y)=1] ans=i=1nj=1mxiyj[gcd(x,y)=1]
根据莫比乌斯反演,得到
a n s = ∑ i = 1 n ∑ j = 1 m ∑ x ∣ i ∑ y ∣ j ∑ d ∣ g c d ( x , y ) μ ( d ) ans=\sum^n_{i=1}\sum^m_{j=1}\sum_{x|i}\sum_{y|j}\sum_{d|gcd(x,y)}\mu(d) ans=i=1nj=1mxiyjdgcd(x,y)μ(d)
直接枚举 d d d,得到
a n s = ∑ i = 1 n ∑ j = 1 m ∑ x ∣ i ∑ y ∣ j ∑ d = 1 m i n ( n , m ) μ ( d ) ∗ [ d ∣ g c d ( x , y ) ] ans=\sum^n_{i=1}\sum^m_{j=1}\sum_{x|i}\sum_{y|j}\sum_{d=1}^{min(n,m)}\mu(d)*[d|gcd(x,y)] ans=i=1nj=1mxiyjd=1min(n,m)μ(d)[dgcd(x,y)]
把只有关 d d d的项移出,得到
a n s = ∑ d = 1 m i n ( n , m ) μ ( d ) ∑ i = 1 n ∑ j = 1 m ∑ x ∣ i ∑ y ∣ j [ d ∣ g c d ( x , y ) ] ans=\sum_{d=1}^{min(n,m)}\mu(d)\sum^n_{i=1}\sum^m_{j=1}\sum_{x|i}\sum_{y|j}[d|gcd(x,y)] ans=d=1min(n,m)μ(d)i=1nj=1mxiyj[dgcd(x,y)]
换一种方式,得到
a n s = ∑ d = 1 m i n ( n , m ) μ ( d ) ∑ x = 1 n ∑ y = 1 m [ d ∣ g c d ( x , y ) ] ⌊ n x ⌋ ⌊ m y ⌋ ans=\sum_{d=1}^{min(n,m)}\mu(d)\sum_{x=1}^n\sum_{y=1}^m[d|gcd(x,y)]\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor ans=d=1min(n,m)μ(d)x=1ny=1m[dgcd(x,y)]xnym
枚举 d x , d y dx,dy dx,dy,得到
a n s = ∑ d = 1 m i n ( n , m ) μ ( d ) ∑ x = 1 ⌊ n d ⌋ ∑ y = 1 ⌊ m d ⌋ ⌊ n d x ⌋ ⌊ m d y ⌋ ans=\sum_{d=1}^{min(n,m)}\mu(d)\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\lfloor\frac{n}{dx}\rfloor\lfloor\frac{m}{dy}\rfloor ans=d=1min(n,m)μ(d)x=1dny=1dmdxndym
再移项,得到
a n s = ∑ d = 1 m i n ( n , m ) μ ( d ) ( ∑ x = 1 ⌊ n d ⌋ ⌊ n d x ⌋ ) ( ∑ y = 1 ⌊ m d ⌋ ⌊ m d y ⌋ ) ans=\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) ans=d=1min(n,m)μ(d)(x=1dndxn)(y=1dmdym)
发现可以用整除分块和线性筛完成,于是


代码

#include <cstdio>
#include <cctype>
#define rr register
using namespace std;
const int N=50000;
int mu[N+1],v[N+1],prime[N+1],cnt,sum[N+1];
inline signed iut(){
    rr int ans=0; rr char c=getchar();
    while (!isdigit(c)) c=getchar();
    while (isdigit(c)) ans=(ans<<3)+(ans<<1)+c-48,c=getchar();
    return ans;
}
inline void prepa(){
    mu[1]=sum[1]=1; 
    for (rr int i=2;i<=N;++i){//线性筛莫比乌斯函数
        if (!v[i]) mu[i]=-1,v[i]=prime[++cnt]=i;
        for (rr int j=1;j<=cnt&&prime[j]*i<=N;++j){
            v[i*prime[j]]=prime[j];
            if (i%prime[j]) mu[i*prime[j]]=-mu[i];
                else break;
        }
    }
    for (rr int i=2;i<=N;++i){
        mu[i]+=mu[i-1];
        for (rr int l=1,r;l<=i;l=r+1)//整除分块
            sum[i]+=((r=(i/(i/l)))-l+1)*(i/l);
    }
}
inline signed min(int a,int b){return (a<b)?a:b;}
signed main(){
    prepa();
    rr int t=iut();
    while (t--){
        rr int n=iut(),m=iut();
        rr int minx=min(n,m);
        rr long long ans=0;
        for (rr int l=1,r;l<=minx;l=r+1){
            r=min(n/(n/l),m/(m/l));
            ans+=(mu[r]-mu[l-1])*1ll*sum[n/l]*sum[m/l];//按照刚刚的式子求出答案
        }
        printf("%lld\n",ans);
    }
    return 0;
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

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

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值