Luogu12522 限りなく灰色へ

~~太好了是无限灰我的大小键跳拍没救了~~

题面

简化题意

给定 $n$ 个点 $(x_i,y_i)$ ，求

\max_{(x,y)}^{[1,X]\times[1,Y]}\sum_{i=1}^n [\gcd(x-x_i,y-y_i)=1]

其中， $\gcd(0,0)=\gcd(0,1)=1,\gcd(0,x)=0\,(x\neq0,x\neq 1)$ .

$n,X,Y\leq 2\times10^3$ . 后文中，为了方便，取值域上界为 $C$ .

题解

60pts

一个能拿满部分分的做法，莫反一步到位：

\begin{aligned} & \max_{(x,y)}^{[1,X]\times[1,Y]}\sum_{i=1}^n [\gcd(x-x_i,y-y_i)=1] \\ =& \max_{(x,y)}^{[1,X]\times[1,Y]}\sum_{i=1}^n \sum_{d\mid\gcd(x-x_i,y-y_i)} \mu(d) \\ =& \max_{(x,y)}^{[1,X]\times[1,Y]}\sum_{d=1}^C \mu(d) \sum_{i=1}^n [d\mid x-x_i][d\mid y-y_i] \end{aligned}

考虑 $d\mid x-x_i \Leftrightarrow x\equiv x_i\pmod d$ （ $y$ 同理）。考虑根号分治，设置阈值 $B$ ：

对于小 $d$ ，更新余数点对 $(x_i\bmod d, y_i\bmod d)$ .
对于大 $d$ ，暴力更新所有符合条件的点对 $(x,y)$ .

那么一个点 $(x,y)$ 的答案就是大 $d$ 的答案加上小 $d$ 对应的余数点对 $(x\bmod d, y\bmod d)$ 的答案。

得到一种实现为：

// not0 为 mu(d)!=0 的 d
for(auto d:not0){
    if(d<=B){
        useD.push_back(d);
        for(int i=1;i<=n;++i)
            mp[Tuple(d,a[i].x%d,a[i].y%d)]+=mu[d];
    }
    else{
        for(int i=1;i<=n;++i)
            for(int x=a[i].x%d;x<=X;x+=d)
                for(int y=a[i].y%d;y<=Y;y+=d)
                    ans[x][y]+=mu[d];
    }
}
for(int x=1;x<=X;++x)
    for(int y=1;y<=Y;++y){
        for(auto d:useD)
            ans[x][y]+=mp[Tuple(d,x%d,y%d)];
        mx=max(mx,ans[x][y]);
    }

关于 $B$ 的设置：对于小 $d$ ，复杂度为 $\mathcal O(Bn\log n + C^2B\log n) = \mathcal O(B(n+C^2)\log n)$ ；对于大 $d$ ，复杂度为 $\mathcal O\left(Bn\cdot\left(\frac{C}{B}\right)^2\right) = \mathcal O\left(\frac{nC^2}{B}\right)$ . 相等时：

B = \sqrt{\frac{nC^2}{(n+C^2)\log n}}

实际取 $B=\min\{C,15\}$ ，60pts 的部分不到 0.25s 就跑出来了。

100pts

上面的做法中我们总是枚举点对，对所有点直接算答案，现在考虑枚举 $x$ ，动态维护 $y$ 的答案。

重新考虑前面莫反的式子，一种做法是：枚举 $x$ ，对于每一个点 $(x_i,y_i)$ ，枚举 $d\mid x-x_i$ ，更新满足 $d\mid y-y_i$ （即 $y\equiv y_i\pmod d$ ）的 $y$ 的答案。沿用 60pts 的根号分支思想，设置阈值 $B$ ：

对于小 $d$ ，设 $cnt_{d,r}$ 表示满足 $d\mid x-x_i$ 且 $y_i\equiv r\pmod d$ 的点 $(x_i,y_i)$ 的个数，用 $cnt_{d,r}\cdot\mu(d)$ 批量更新。
对于大 $d$ ，暴力更新。

需要注意 $x=x_i$ 的情况细节。得到一种实现为：

for(int x=1;x<=X;++x){
    fill(ans+1,ans+Y+1,count_if(a+1,a+n+1,[&](Point p){return p.x!=x;}));
    for(int i=1;i<=n;++i) if(x==a[i].x)
        ++ans[a[i].y],++ans[a[i].y-1],++ans[a[i].y+1];
    // x=x_i 的细节
    for(int i=1;i<=n;++i) if(x!=a[i].x){
        fact(abs(x-a[i].x));
        // 此函数求得 not0fac，即满足 mu(d)!=0 且 d>1 的 deltaX 的约数 d
        for(auto d:not0fac){
            int r=a[i].y%d;
            if(d<=B)
                cnt[d][r]+=mu[d];
            else{
                for(int y=r;y<=Y;y+=d)
                    ans[y]+=mu[d];
            }
        }
    }
    for(auto d:not0){
        if(d>B) break;
        for(int r=0;r<d;++r) if(cnt[d][r]){
            for(int y=r;y<=Y;y+=d)
                ans[y]+=cnt[d][r];
            cnt[d][r]=0;
        }
    }
    mx=max(mx,*max_element(ans+1,ans+Y+1));
}

关于 $B$ 的设置，设 $m_d$ 为 not0fac 的大小：对于小 $d$ ，复杂度为 $\mathcal O(Cnm_d+C^2B)$ ；对于大 $d$ ，复杂度为 $\mathcal O\left(\frac{C^2nm_d}{B}\right)$ . 相等时：

B = \sqrt{nm_d}

注意到 $2\times3\times5\times7<2000<2\times3\times5\times7\times11$ ，故 $m_d$ 最大为 $2^4=16$ . 实际取 $B=\min\{C,180\}$ 可过。

~~实际上我代码太丑了，怎么取都跑巨慢。~~

代码

#include<bits/stdc++.h>
// #define int long long
#ifdef DEBUG
#define debug(...) fprintf(stderr,__VA_ARGS__)
#else
#define debug
#endif
using namespace std;
bool Begin;
const int max_n=2003;
int read(){
    int x=0;bool w=0;char c=getchar();
    while(!isdigit(c)) w|=c=='-',c=getchar();
    while(isdigit(c)) x=x*10+(c^48),c=getchar();
    return w?-x:x;
}

int n,X,Y,C;

vector<int> P;
bool isp[max_n];
int mu[max_n];
void sieve(int n){
    memset(isp,1,sizeof(isp));
    isp[1]=0;
    mu[1]=1;
    for(int i=2;i<=n;++i){
        if(isp[i]){
            P.push_back(i);
            mu[i]=-1;
        }
        for(auto p:P){
            if(i*p>n) break;
            isp[i*p]=0;
            if(i%p==0){
                mu[i*p]=0;
                break;
            }
            mu[i*p]=-mu[i];
        }
    }
}

struct Point{
    int x,y;
    Point(int x=0,int y=0):x(x),y(y){}
}a[max_n];

vector<int> not0;

int ans[max_n],cnt[max_n][max_n];

vector<int> not0fac;
void fact(int x){
    vector<int> facP;
    for(auto p:P){
        if(x%p==0){
            facP.push_back(p);
            while(x%p==0) x/=p;
            if(x==1) break;
        }
    }
    not0fac.clear();
    int n=facP.size(),S=1<<n;
    for(int i=1;i<S;++i){
        int d=1;
        for(int j=0;j<n;++j)
            if(i>>j&1) d*=facP[j];
        not0fac.push_back(d);
    }
}

bool End;
#define File ""
signed main(){
    #ifndef ONLINE_JUDGE
    freopen(File ".in","r",stdin),
    freopen(File ".out","w",stdout);
    #endif
    debug("Memory: %lf\n",(&Begin-&End)/1024.0/1024);
    X=read(),Y=read(),n=read(),C=max(X,Y);
    for(int i=1;i<=n;++i){
        int x=read(),y=read();
        C=max({C,x,y});
        a[i]=Point(x,y);
    }
    int B=min(C,180),mx=0;
    sieve(C);
    for(int d=2;d<=C;++d) if(mu[d])
        not0.push_back(d);
    for(int x=1;x<=X;++x){
        fill(ans+1,ans+Y+1,count_if(a+1,a+n+1,[&](Point p){return p.x!=x;}));
        for(int i=1;i<=n;++i) if(x==a[i].x)
            ++ans[a[i].y],++ans[a[i].y-1],++ans[a[i].y+1];
        for(int i=1;i<=n;++i) if(x!=a[i].x){
            fact(abs(x-a[i].x));
            for(auto d:not0fac){
                int r=a[i].y%d;
                if(d<=B)
                    cnt[d][r]+=mu[d];
                else{
                    for(int y=r;y<=Y;y+=d)
                        ans[y]+=mu[d];
                }
            }
        }
        for(auto d:not0){
            if(d>B) break;
            for(int r=0;r<d;++r) if(cnt[d][r]){
                for(int y=r;y<=Y;y+=d)
                    ans[y]+=cnt[d][r];
                cnt[d][r]=0;
            }
        }
        mx=max(mx,*max_element(ans+1,ans+Y+1));
    }
    printf("%d\n",mx);
    return 0;
}