【题解】「Ynoi2018」GOSICK [*hard]

这题为什么是 hard？显然因为这题贼卡常 = =

其实要处理的就是前缀的约数个数和和倍数个数和（后缀同样）。

等于要做的就是对于 (a_i) ，询问 (a_{1,2,cdots,i-1}) 中有多少是 (i) 的约数，又有多少是 (i) 的倍数。（对于 (a_i) 相等的不计算）。

分别考虑约数和倍数怎么做。

• 约数的做法

因为一个数的约数最多就那么多个 ......

所以直接暴力遍历约数就好了

• 倍数的做法

这会儿不能暴力遍历倍数了，因为如果给你一堆的 1 的话会直接暴毙 = =

直接查询，加入一个数的时候遍历其约数即可。

• 做一遍莫队，计算前后缀和贡献

这个直接做就好了吧？没什么复杂度的问题。

• 修正答案

注意到这里按照最开始的办法做倍数是可以的，因为查询是 (O(1)) 的。

但是约数就不行。考虑怎么将约数的查询弄成 (O(1)) 或者找更好的办法。

（然后就不是很会做了 = =

考虑奇怪的做法：

• (a_i> L) ：暴力跳倍数。
• (a_ileq L)：直接用二进制状态记录。

查询的时候，先处理每个数的 (leq L) 的所有因子组成的状态，然后直接对标查询。

（ps：这个做法常数太大，过不去

#define debug

const int N=5e5+5;
const int M=7e6+5;

int n,m,a[N];
ll ans[N];
short facsta[N][4];
struct Query {int l,r,id;} q[N];

int edge_cnt,head[N<<2];
struct Edge {int nxt,to;} G[M];
inline void addedge(int u,int v) {G[++edge_cnt]=(Edge){head[u],v},head[u]=edge_cnt;}

namespace InitFactor { // {{{ Init Factor
    const int L=N-5;

    int cnt,d[N],prime[N];
    bool vis[N];
    pii hep[N];

    inline void sieve() { // Get d
        lep(i,2,L) {
            if(!vis[i]) prime[++cnt]=i,d[i]=i;
            for(int j=i+i;j<=L;j+=i)
                if(!vis[j]) d[j]=i,vis[j]=true;
        }
    }
    void dfs(int step,int rt,int res) {
        if(step>cnt) return addedge(rt,res),void();
        lep(i,0,hep[step].se) dfs(step+1,rt,res),res*=hep[step].fi;
    }
    inline void solve() {
        sieve();
        lep(i,1,L) {
            int now=i,las=-1,tot=0; cnt=0;
            while(now!=1) {
                if(d[now]==las) ++tot;
                else {if(~las) hep[++cnt]=mkp(las,tot); las=d[now],tot=1;}
                now/=d[now];
            }
            hep[++cnt]=mkp(las,tot),dfs(1,i,1);
        }
        lep(i,1,L) lep(j,0,7) lep(k,0,3) facsta[i][k]|=(i%(j+k*8+1)==0)<<j;
    }
} // }}}

namespace InitPreSuf { // {{{ Get Pre1,Suf1,Pre2,Suf2
    ll pre1[N],suf1[N],pre2[N],suf2[N];
    int sta[N];

    inline void solve() {
        CLEAR(sta); lep(i,1,n) { // pre1
            for(int e=head[a[i]];e;e=G[e].nxt) pre1[i]+=sta[G[e].to];
            ++sta[a[i]],pre1[i]+=pre1[i-1];
        }
        CLEAR(sta); rep(i,n,1) { // suf1
            for(int e=head[a[i]];e;e=G[e].nxt) suf1[i]+=sta[G[e].to];
            ++sta[a[i]],suf1[i]+=suf1[i+1];
        }
        CLEAR(sta); lep(i,1,n) { // pre2
            pre2[i]=pre2[i-1]+sta[a[i]];
            for(int e=head[a[i]];e;e=G[e].nxt) ++sta[G[e].to];
        }
        CLEAR(sta); rep(i,n,1) { // suf2
            suf2[i]=suf2[i+1]+sta[a[i]];
            for(int e=head[a[i]];e;e=G[e].nxt) ++sta[G[e].to];
        }
    }
}
using InitPreSuf::pre1;
using InitPreSuf::pre2;
using InitPreSuf::suf1;
using InitPreSuf::suf2;
// }}}

namespace FixAnswer { // {{{ Fix Answer
    struct Node {int i,l,r,id,typ;} sta1[N<<2],sta2[N<<2];
    int cnt1,cnt2;

    // {{{ Struct
    int sum[4][256],sta[N];

    inline void update(int x) {
        if(x<=32) {
            int id=(x-1)/8; x=1<<(x-id*8-1);
            lep(S,0,255) if(S&x) ++sum[id][S];
        } else for(int i=x;i<N;i+=x) ++sta[i];
    }
    inline int query(int x) {
        return sum[0][facsta[x][0]]+sum[1][facsta[x][1]]+
               sum[2][facsta[x][2]]+sum[3][facsta[x][3]]+sta[x];
    }
    // }}}

    inline void solve() {
        std::sort(sta1+1,sta1+1+cnt1,[](Node x,Node y){return x.i<y.i;});
        std::sort(sta2+1,sta2+1+cnt2,[](Node x,Node y){return x.i>y.i;});
        int l;
        
        l=1,CLEAR(sum),CLEAR(sta); lep(i,1,n) {
            update(a[i]);
            while(l<=cnt1&&sta1[l].i==i) {
                lep(j,sta1[l].l,sta1[l].r) ans[sta1[l].id]+=query(a[j])*sta1[l].typ;
                ++l;
            }
        }
        l=1,CLEAR(sum),CLEAR(sta); rep(i,n,1) {
            update(a[i]);
            while(l<=cnt2&&sta2[l].i==i) {
                lep(j,sta2[l].l,sta2[l].r) ans[sta2[l].id]+=query(a[j])*sta2[l].typ;
                ++l;
            }
        }
        
        l=1,CLEAR(sta); lep(i,1,n) {
            for(int e=head[a[i]];e;e=G[e].nxt) ++sta[G[e].to];
            while(l<=cnt1&&sta1[l].i==i) {
                lep(j,sta1[l].l,sta1[l].r) ans[sta1[l].id]+=sta[a[j]]*sta1[l].typ;
                ++l;
            }
        }
        l=1,CLEAR(sta); rep(i,n,1) {
            for(int e=head[a[i]];e;e=G[e].nxt) ++sta[G[e].to];
            while(l<=cnt2&&sta2[l].i==i) {
                lep(j,sta2[l].l,sta2[l].r) ans[sta2[l].id]+=sta[a[j]]*sta2[l].typ;
                ++l;
            }
        }

        lep(i,1,m) ans[i]+=ans[i-1];
    }
}
using FixAnswer::Node;
using FixAnswer::sta1;
using FixAnswer::sta2;
using FixAnswer::cnt1;
using FixAnswer::cnt2;
// }}}

namespace GetAnswer { // {{{ Get Answer
    inline void solve() {    
        int l=1,r=0;
        lep(i,1,m) {
            if(r<q[i].r) {
                ans[i]+=pre1[q[i].r]-pre1[r]+pre2[q[i].r]-pre2[r];
                if(l>1) sta1[++cnt1]=(Node){l-1,r+1,q[i].r,i,-1}; r=q[i].r;
            }
            if(r>q[i].r) {
                ans[i]-=pre1[r]-pre1[q[i].r]+pre2[r]-pre2[q[i].r];
                if(l>1) sta1[++cnt1]=(Node){l-1,q[i].r+1,r,i,1}; r=q[i].r;
            }
            if(l<q[i].l) {
                ans[i]-=suf1[l]-suf1[q[i].l]+suf2[l]-suf2[q[i].l];
                if(r<n) sta2[++cnt2]=(Node){r+1,l,q[i].l-1,i,1}; l=q[i].l;
            }
            if(l>q[i].l) {
                ans[i]+=suf1[q[i].l]-suf1[l]+suf2[q[i].l]-suf2[l];
                if(r<n) sta2[++cnt2]=(Node){r+1,q[i].l,l-1,i,-1}; l=q[i].l;
            }
        }
    }
} // }}}

int pos[N];
int main() {
    InitFactor::solve();

    IN(n,m);
    lep(i,1,n) IN(a[i]);
    lep(i,1,m) IN(q[i].l,q[i].r),q[i].id=i;

    int B=n/sqrt(m*2/3+1);
    std::sort(q+1,q+1+m,[&](Query a,Query b){return (a.l/B)^(b.l/B)?a.l<b.l:(((a.l/B)&1)?a.r<b.r:a.r>b.r);});
    lep(i,1,m) pos[q[i].id]=i;

    InitPreSuf::solve();
    GetAnswer::solve();
    FixAnswer::solve();

    lep(i,1,m) printf("%lld
",ans[pos[i]]+q[pos[i]].r-q[pos[i]].l+1);
    return 0;
}

这个想法固然十分巧妙，但是过不去。

同时这个做法又很长 = =

考虑一种更简单的做法：根号分治：

• 对于 (a_i>sqrt{n}) ：暴力跳倍数。
• 对于 (a_ileq sqrt{n}) ：预处理 (sum_{x,i}) 表示 (a_{1,2,cdots,i}) 中 (x) 的倍数个数，查询的话就遍历 ([1,sqrt{n}]) 中的每个数，然后看看在这段区间中的倍数个数即可。

需要注意的是虽然莫队的移动量是 (O(nsqrt{m})) ，但是这个是对一段区间查询，因此查询总数是 (O(m)) 的，所以说虽然一次查询的复杂度是 (O(sqrt{n}))，这一块复杂度也还是 (O(msqrt{n})) 的。

然后没过几天 lxl 就改空间限制了 = = ，根本开不下。

首先有一种空间 (O(n)) 的做法：遍历 ([1,sqrt{n}]) 中的每一个数，然后做前缀和，接着在遍历一遍所有的询问，处理贡献。不难发现这个做法常数比较大。

接着优化：将数字分成 (8) 个一组，这样空间还是 (O(n)) 的（开 (8) 个 (sum) 数组），然后将这 (8) 个数连着处理即可，这样就优秀不少。

还需要注意的是，因为常熟问题，暴力跳倍数的常数很小，因此可以做一些平衡：我实现的时候是将这个阈值设为了 (80) 。

（除了 FixAnswer，其他部分和上面一模一样 = =

const int T=80;

namespace FixAnswer { // {{{ Fix Answer
    struct Node {int i,l,r,id,typ;} seq1[N],seq2[N];
    int L,cnt1,cnt2,sta[N],tot[T+2],sum[8][N];

    inline void update(int x) {
        for(int e=head[x];e;e=G[e].nxt) ++sta[G[e].to];
        if(x>T) for(int i=x;i<N;i+=x) ++sta[i];
    }
    inline void query(Node now) {
        lep(i,now.l,now.r) ans[now.id]+=sta[a[i]]*now.typ;
    }
    inline void calc(Node now) {
        int res=0;
        lep(_,0,7) res+=tot[L+_]*(sum[_][now.r]-sum[_][now.l-1]);
        ans[now.id]+=res*now.typ;
    }

    inline void solve() {
        std::sort(seq1+1,seq1+1+cnt1,[](Node x,Node y){return x.i<y.i;});
        std::sort(seq2+1,seq2+1+cnt2,[](Node x,Node y){return x.i>y.i;});
        int l;

        l=1,CLEAR(sta); lep(i,1,n) {update(a[i]); while(l<=cnt1&&seq1[l].i==i) query(seq1[l]),++l;}
        l=1,CLEAR(sta); rep(i,n,1) {update(a[i]); while(l<=cnt2&&seq2[l].i==i) query(seq2[l]),++l;}
        
        lep(_,1,9) {
            L=(_-1)*8+1;
            lep(j,0,7) {
                const int num=L+j;
                lep(i,1,n) sum[j][i]=sum[j][i-1]+(a[i]%num==0);
            }

            l=1,CLEAR(tot); lep(i,1,n) {if(a[i]<=T) ++tot[a[i]]; while(l<=cnt1&&seq1[l].i==i) calc(seq1[l]),++l;}
            l=1,CLEAR(tot); rep(i,n,1) {if(a[i]<=T) ++tot[a[i]]; while(l<=cnt2&&seq2[l].i==i) calc(seq2[l]),++l;}
        }
        
        lep(i,1,m) ans[i]+=ans[i-1];
    }
} using FixAnswer::Node,FixAnswer::seq1,FixAnswer::seq2,FixAnswer::cnt1,FixAnswer::cnt2; // }}}