题解 P5624 [Celeste-A]Black Moonrise 暗夜月升

· · 题解

update: 修正了错误代码。

题目传送门

前置知识：

欧拉函数，分块。

题意：

分析：

&\sum\limits_{i=l}^r\sum\limits_{j=l}^r\gcd(a_i,a_j)\\ =&\sum\limits_{i=l}^r\sum\limits_{j=l}^r\sum\limits_{d\mid a_i\land d\mid a_j}\varphi(d)&(1)\\ =&\sum\limits_{d=1}^{\max(a)}\varphi(d)\sum\limits_{i=l}^r\sum\limits_{j=l}^r[d\mid a_i\land d\mid a_j]\\ =&\sum\limits_{d=1}^{\max(a)}\varphi(d)(\sum\limits_{i=l}^r[d\mid a_i])^2 \end{aligned}

其中 (1) 用到了欧拉函数的性质 n=\sum\limits_{d\mid n}\varphi(d)（也可用莫反推出一个新函数 f(d)，实际上 f 就是 \varphi）

官方解法中使用了莫队维护 \sum\limits_{i=l}^r[d\mid a_i]（不会莫队可以跳过这一段）：由于 a_i 随机生成，故其因数个数期望较少（另一篇题解中认为是一个常数 \alpha，本题解认为是值域内的总因子个数 v\ln v 除以值域 v，即 \ln v），移动莫队指针时可以直接枚举所有 a_i 的因数 d 并更新对应的 \sum\limits_{i=l}^r[d\mid a_i]。期望复杂度 \mathcal O(n\sqrt n\ln v)，惊人地可过。

但翻了 AC 代码，发现本题实际上有一种期望 \mathcal O((n+q)\sqrt n) 的分块写法。

首先对值域分块，所有不大于 V 的 d 直接暴力求出 [d\mid a_i] 的前缀和，询问时 \mathcal O(V) 即可处理这一部分的数，总复杂度 \mathcal O((n+q)V)。

对于大于 V 的 d，由于值比较大，所以随机情况下 a 中 d 的倍数较少。假设所有为 d 的倍数的 a 为 a_{k_1},a_{k_2},\dots,a_{k_m}。当询问左端点 l\in(k_{x-1},k_x]，右端点 r\in[k_y,k_{y+1}) 时，\sum\limits_{i=l}^r[d\mid a_i] 即为 y-x+1；对答案贡献为 \varphi(d)(y-x+1)^2。对于每一个询问，显然不能循环每一个 d 计算 x,y。但注意到因为 m 较小，有序对 (x,y) 只有 \mathcal O(m^2) 种情况，可以从这里入手。

不妨将询问按右端点排序，对于 y 相同的询问，d 的贡献只与左端点有关。而左端点落在一个区间内时，贡献都是一样的。考虑差分，则 d 只对 y+1 个位置有影响。每次 y 改变只需修改这几处的值即可。总修改次数为 \mathcal O(m^2)。

可以证明 m^2 的期望为 \mathcal O(\dfrac{n^2}{d^2})（证明 link）。则对于所有 d，总修改次数期望为 \mathcal O(\sum\limits_{d=V}^{\max(a)}\dfrac{n^2}{d^2})=\mathcal O(\dfrac{n^2}V)。

对每个询问需要计算出差分数组的前缀和。考虑分块，使修改 \mathcal O(1)，查询 \mathcal O(\sqrt n)。总时间复杂度 \mathcal O((n+q)V+\dfrac{n^2}V+q\sqrt n)，取 V=\sqrt n，复杂度即为 \mathcal O((n+q)\sqrt n)。

思路：

  1. 求欧拉函数。

  2. 对值域分块，分别处理。

具体说明值域分块后大于 V 的 d 如何处理。

开 \max(a) 个 vector，存下每个 d 的 k_i。考虑右端点所在区间由 [k_{y_0},k_{y_0+1}) 到 [k_{y_0+1},k_{y_0+2})，差分数组如何变化：\forall x\in[1,y_0+1]，左端点在 (k_{x-1},k_x] 的贡献由 \varphi(d)(y_0-x+1)^2 变为 \varphi(d)(y_0-x+2)^2，增加 \varphi(d)(2(y_0-x+1)+1)。

在差分数组上，体现为位置 1 增加 \varphi(d)(2y_0+1)，位置 k_x+1\ (x\in[1,y_0]) 减少 2\varphi(d)，位置 k_{y_0}+1 减少 \varphi(d)。

分块求和，记得处理对应位置块和。

#include <bits/stdc++.h>
#define rep(i, l, r) for(int i=l, _=r; i<=_; ++i)
typedef long long ll;
using namespace std;
inline int read() {
    int res=0; bool k=1; char ch;
    while(!isdigit(ch=getchar())) k=ch^45;
    do res=(res<<1)+(res<<3)+(ch&15); while(isdigit(ch=getchar()));
    return k? res: -res;
}
const int mN=1e5+9, mP=1e4, mL=320;
int n, q, L, V, a[mN], sum[mN][mL];
ll ans[mN];

int pri[mP], cntp, phi[mN];
inline void init(int n) {
    phi[1]=1;
    for(int i=2, j, u; i<=n; ++i) {
        if(!phi[i]) pri[++cntp]=i, phi[i]=i-1;
        for(j=1; (u=pri[j]*i)<=n && i/pri[j]*pri[j]^i; ++j) phi[u]=phi[i]*(pri[j]-1);
        if(u<=n) phi[u]=phi[i]*pri[j];
    }
}

struct Que {int x, y, id;} qn[mN];
bool operator < (Que q1, Que q2) {return q1.y<q2.y;}

ll b[mN], sumb[mL]; //b 为差分数组，sumb 为分块的差分数组
vector<int> num[mN];    //用于存 ki
inline void sol(int x, int p) {
    if(x<=V) return ++sum[p][x], void();    //小于 V，直接前缀和
    int y0=(int) num[x].size(), t=phi[x];
    b[1]+=(ll) (y0<<1|1)*t, sumb[0]+=(ll) (y0<<1|1)*t;
    b[p+1]-=t, sumb[(p+1)/L]-=t;
    t<<=1;
    for(int i=y0-1; i>=0; --i) b[num[x][i]+1]-=t, sumb[(num[x][i]+1)/L]-=t;
    num[x].push_back(p);
}

inline ll sqr_(int x) {return (ll) x*x;}

int main() {
    init(1e5), n=read(), q=read(), L=sqrt(n), V=sqrt(1e5);
    if(L<2) L=2;
    rep(i, 1, n) a[i]=read();
    rep(i, 1, q) qn[i].x=read(), qn[i].y=read(), qn[i].id=i;
    sort(qn+1, qn+q+1);

    rep(i, 1, q) {
        rep(r, qn[i-1].y+1, qn[i].y) {  //扩展 y
            memcpy(sum[r], sum[r-1], sizeof sum[r]);
            for(int j=1; j*j<=a[r]; ++j) if(a[r]%j==0) {
                sol(j, r);
                if(j*j^a[r]) sol(a[r]/j, r);
            }
        }
        int x=qn[i].x, y=qn[i].y, id=qn[i].id;
        rep(j, 0, x/L-1) ans[id]+=sumb[j];  //大于 V，值域分块
        rep(j, x/L*L, x) ans[id]+=b[j];
        rep(v, 1, V) ans[id]+=phi[v]*sqr_(sum[y][v]-sum[x-1][v]);   //小于 V，直接前缀和差分
    }
    rep(i, 1, q) printf("%lld\n", ans[i]);
    return 0;
}