更优的区间互质对个数算法

目前我可以做到 $O(n^{1.5}\log n)$ 时间， $O(n)$ 空间。

题目描述：

给你一个长度为 $n$ 的排列 $a$ ， $m$ 次询问，每次询问区间 $[l,r]$ 内互质对个数。

互质对定义： 一个有序二元组 $(i,j)$ 为互质对，当且仅当 $i<j$ 且 $\gcd(a_i,a_j)=1$ 。

区间内二元组定义： 一个有序二元组 $(i,j)$ 在区间 $[l,r]$ 内，当且仅当 $l \le i<j \le r$ 。

$n,m \le 5\times 10^4$ ，

时间限制 $1s$ ，空间限制 $512MB$ 。

首先使用莫队暴力移动需要每次得到一个区间内有多少个数字于已知数字互质。这个可以预处理质数和因数然后 $2^\omega(n)\omega(n)$ 容斥，其中 $\omega(5 \times 10^4)=7$ 。总时间复杂度 $n^{1.5}2^{\omega(n)}\omega(n)$ ，即 $O(896n^{1.5})$ 。

容斥这条路已经走不通了。我们考虑莫队二次离线，那么问题可以转化成往一个集合中加入 $x$ ，以及查询集合中有多少个数与 $x$ 互质。

也就是求 $\displaystyle\sum_{i}[\gcd(i,x)=1]$ 。

莫比乌斯反演之后可以得到 $\displaystyle\sum_{d} \mu(d)\sum_{d\mid i} 1$ 。

令 $cnt_d$ 为集合中有多少个数 $i$ 满足 $d\mid i$ 。

那么我们就可以做到往集合里加入一个数的时间和查询的时间均为 $O(d(n))$ ，其中 $d(n)$ 为约数个数，在 $n=5\times 10^4$ 时与 $2^{\omega(n)}$ 相当。总时间复杂度为 $O(n^{1.5}d(n))$ ，即 $128n^{1.5}$ 。

由于莫队二次离线有 $O(n)$ 次插入， $O(n^{1.5})$ 次查询，并且插入是每个数字插入一次，总复杂度 $\displaystyle\sum_{i \le n}d(i)=O(n\log n)$ 。尝试 $O(1)$ 查询。那么就需要维护答案数组 $ans_i$ 。具体的，对于加入一个数字 $x$ 的操作，需要对于每一个 $d\mid x$ ，对所有 $d\mid y$ ， $ans_y :=ans_y+\mu(y)$ 。容易得到这样做插入复杂度是 $O(4\times 10^{10})$ 级别，不可接受。可以发现大多数复杂度都是较小的 $d$ 贡献的，于是考虑根号分治。

设置一个阈值 $B$ ，当 $d\ge B$ 时，暴力更新 $ans$ 数组。当 $d<B$ 时，采用 $cnt$ 数组进行更行。取 $B=16$ ，可得总复杂度 $O(n^{1.5}\log n)$ 。

代码：

#include <bits/stdc++.h>
using namespace std;
const int N=50010,B=230,T=16;
int n,m,a[N],ans[N],pre[N],c[N],d[N],prel[N];
struct node{
    int l,r,id,ans;
}q[N];
struct Node{
    int l,r,opt,id;
}; vector<Node>v[N];
bool cmp(node x,node y){
    return x.l/B==y.l/B?x.r/B<y.r/B:x.l/B<y.l/B;
}
namespace MU{
    int book[N],p[N],cnt,mu[N];
    vector<int>pl[N];
    void get_prime(int x){
        mu[1]=1;
        for(int i=2;i<=x;i++){
            if(!book[i])p[++cnt]=i,mu[i]=-1;
            for(int j=1;j<=cnt&&i*p[j]<=x;j++){
                book[i*p[j]]=1;
                if(i%p[j]==0)break;
                mu[i*p[j]]=-mu[i];
            }
        }
        for(int i=1;i<=x;i++){
            for(int j=i;j<=n;j+=i){
                pl[j].push_back(i);
            }
        }
    }
}using namespace MU;
int get(int x){
    int ans=d[x];
    for(int u:pl[x]){
        if(u>=T)break;
        ans+=c[u];
    }
    return ans;
}
int main(){
    cin>>n>>m;
    get_prime(n);
    for(int i=1;i<=n;i++)cin>>a[i];
    for(int i=1;i<=m;i++){
        cin>>q[i].l>>q[i].r;
        q[i].id=i;
    }
    sort(q+1,q+1+m,cmp);
    int l=1,r=0;
    for(int i=2;i<=n;i++){
        for(int u:pl[a[i-1]]){
            c[u]+=mu[u];
        }
        for(int u:pl[a[i]]){
            pre[i]+=c[u];
        }
        prel[i]=pre[i];
    }
    for(int i=1;i<=n;i++)if(a[i]==1)prel[i]++;
    memset(c,0,sizeof(c));
    for(int i=1;i<=m;i++){
        if(l<q[i].l)v[r].push_back({l,q[i].l-1,-1,i});
        while(l<q[i].l)q[i].ans+=prel[l],l++;
        if(l>q[i].l)v[r].push_back({q[i].l,l-1,1,i});
        while(l>q[i].l)l--,q[i].ans-=prel[l];
        if(r<q[i].r)v[l-1].push_back({r+1,q[i].r,-1,i});
        while(r<q[i].r)r++,q[i].ans+=pre[r];
        if(r>q[i].r)v[l-1].push_back({q[i].r+1,r,1,i});
        while(r>q[i].r)q[i].ans-=pre[r],r--;
    }
    for(int i=1;i<=n;i++){
        for(int u:pl[a[i]]){
            if(u>=T){
                for(int j=u;j<=n;j+=u)d[j]+=mu[u];
            }
            else{
                c[u]+=mu[u];
            }
        }
        for(Node j:v[i]){
            for(int k=j.l;k<=j.r;k++){
                q[j.id].ans+=j.opt*get(a[k]);
            }
        }
    }
    for(int i=1;i<=m;i++)q[i].ans+=q[i-1].ans;
    for(int i=1;i<=m;i++)ans[q[i].id]=q[i].ans;
    for(int i=1;i<=m;i++)cout<<ans[i]<<"\n";
    return 0;
}