[POI2007]ZAP-Queries

Posted 2021-10-31 a1b3c7d9

[POI2007]ZAP-Queries

求n组\(\sum_{i=1}^a\sum_{j=1}^b(gcd(i,j)==d)\), 1≤n≤50 000,1≤d≤a,b≤50 000。

解

不难看出是约数组合计数问题，而解决该问题常用思路有容斥原理和Mobius反演。

法一：容斥原理

尽可能特殊化问题，令\(a/=d,b/=d\)，原问题等价于在这个范围里寻找互质的数的对数，而容斥原理关键在于构造不合理情况，即有别的约数，于是设\(D(a,b,i)\)表示a，b范围里有公约数i的个数，不难得知\(D(a,b,i)=[a/i][b/i]\),于是由容斥原理有
\[ans=\sum_{i=1}^{min(a,b)}\mu(i)D(a,b,i)=\sum_{i=1}^{min(a,b)}\mu(i)[a/i][b/i]\]
预处理\(\mu\)函数前缀和，进行整出分块即可。

法二：Mobius反演

设
\[f(d)=\sum_{i=1}^a\sum_{j=1}^b(gcd(i,j)==d)\]
\[F(d)=\sum_{i=1}^a\sum_{j=1}^b(d|gcd(i,j))=[a/d][b/d]\]
由Mobius反演定理我们有
\[f(d)=\sum_{d|x}F(x)\mu(x/d)=\sum_{d|x}[a/x][b/x]\mu(x/d)\]
于是
\[ans=f(d)=\sum_{d|x}[a/x][b/x]\mu(x/d)=\]
\[\sum_{x=1}^{min([a/d],[b/d)]}[a/x][b/x]\mu(x/d)\]

对式子进行整除分块即可。

参考代码

#include <iostream>
#include <cstdio>
#define il inline
#define ri register
#define swap(x,y) x^=y^=x^=y
using namespace std;
bool check[50001];
int prime[5137],pt,mb[50001];
il void read(int&);
il int min(int,int);
void pen(int),prepare(int);
int main(){
    int lsy,a,b,n,i,j,ans;read(lsy);
    prepare(50000);while(lsy--){
        read(a),read(b),read(n),ans&=0;
        a/=n,b/=n;if(a>b)swap(a,b);
        for(i=1;i<=a;i=j+1)
            j=min(a/(a/i),b/(b/i)),
                ans+=(a/i)*(b/i)*(mb[j]-mb[i-1]);
        pen(ans),putchar('\n');
    }
    return 0;
}
void prepare(int n){
    int i,j;check[1]|=mb[1]|=true;
    for(i=2;i<=n;++i){
        if(!check[i])mb[i]=-1,prime[++pt]=i;
        for(j=1;j<=pt&&i*prime[j]<=n;++j){
            check[i*prime[j]]|=true;
            if(!(i%prime[j]))break;
            mb[i*prime[j]]=-mb[i];
        }
    }for(i=1;i<=n;++i)mb[i]+=mb[i-1];
}
il int min(int a,int b){
    return a<b?a:b;
}
void pen(int x){
    if(x>9)pen(x/10);putchar(x%10+48);
}
il void read(int &x){
    x&=0;ri char c;while(c=getchar(),c<'0'||c>'9');
    while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+(c^48),c=getchar();
}