[HAOI2011]Problem b

[HAOI2011]Problem b

题目

解法

定义ans(N,M)表示\[\sum_{i=1}^N\sum_{j=1}^M[gcd(i,j)=k]\]
容斥一下答案就是\(ans(b,d)-ans(b,c-1)-ans(a-1,d)+ans(a-1,c-1)\)

我们再考虑单个\(ans(N,M)\)怎么求

先把k提出来
\[ans(N,M)=\sum_{i=1}^\frac{N}{k}\sum_{j=1}^\frac{M}{k}[gcd(i,j)=1]\]
然后考虑一下莫比乌斯反演
定义
\[f(x)=\sum_{i=1}^\frac{N}{k}\sum_{j=1}^\frac{M}{k}[gcd(i,j)=x]\]
\[g(x)=\sum_{x|d}f(d)\]
那么
\[g(x)=\sum_{x|d}\sum_{i=1}^\frac{N}{k}\sum_{j=1}^\frac{M}{k}[gcd(i,j)=d]\]
然后把第一个\(sigma\)去掉
\[g(x)=\sum_{i=1}^\frac{N}{k}\sum_{j=1}^\frac{M}{k}[gcd(i,j)\%x=0]\]
然后把\(x\)提出来去掉\(gcd\)的影响
\[g(x)=\sum_{i=1}^\frac{N}{kx}\sum_{j=1}^\frac{M}{kx}[gcd(i,j)\%1=0]\]
显然
\[g(x)=\frac{N}{kx}\frac{M}{kx}\]
根据莫比乌斯反演定理
\[f(x)=\sum_{x|d}μ(\frac{d}{x})g(d)\]
\[f(x)=\sum_{x|d}μ(\frac{d}{x})\frac{N}{kd}\frac{M}{kd}\]
而\(ans(N,M)=f(1)\)
\[ans(N,M)=\sum_{d=1}^Nμ(d)\frac{N}{kd}\frac{M}{kd}\]
前缀和加数论分块\(O(\sqrt n)\)算出每一次的结果。

然后为了方便统计答案
我们发现当\(kd>N\)时，是没有贡献的
于是乎
\[ans(N,M)=\sum_{d=1}^\frac{N}{k}μ(d)\frac{N}{kd}\frac{M}{kd}\]
所以如果先把\(N,M\)都除掉\(k\)那么
\[ans(N,M)=\sum_{d=1}^Nμ(d)\frac{N}{d}\frac{M}{d}\]
不要忘记容斥,这还不是最终答案。

完整代码

#include<bits/stdc++.h>
#define il inline
#define ll long long
#define rg register
using namespace std;
const int N=1e7;
ll pre[N],prime[N],isprime[N],s;
void work(){
    for(rg int i=2;i<=50000;++i){
        if(!isprime[i])pre[i]=-1,prime[++s]=i;
        for(rg int j=1;j<=s&&i*prime[j]<=50000;++j){
            isprime[i*prime[j]]=1;
            pre[i*prime[j]]=-pre[i];
            if(i%prime[j]==0){
                pre[i*prime[j]]=0;
                break;
            }
        }
    }
    pre[1]=1;
    for(rg int i=1;i<=50000;++i)
            pre[i]+=pre[i-1];
}
il ll js(int x,int y){
    if(x>y)swap(x,y);
    if(x==0)return 0;
    ll ans=0;
    for(rg int l=1,r;l<=x;l=r+1){
        r=min(x/(x/l),y/(y/l));
        ans+=(pre[r]-pre[l-1])*(x/l)*(y/l);
    }
    return ans;
}
int main(){
    work();
    int T;
    cin>>T;
    while(T--){
        int a,b,c,d,k;
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        b/=k,d/=k;
        if(a%k)a=a/k+1;
        else a/=k;
        if(c%k)c=c/k+1;
        else c/=k;
        printf("%lld\n",js(b,d)-js(b,c-1)-js(a-1,d)+js(a-1,c-1));
    }
    return 0;
}

转载于:https://www.cnblogs.com/ljq-despair/p/8683965.html