Luogu P10704 救赎（Redemption）

Luogu P10704 救赎（Redemption）

这个题看起来凉了，所以不再去投题解，在此记录一下 \(O(n+m^{2/3})\) 的做法。
算是 Dirichlet 双曲线法的三维情形的应用示例。

记 \(\displaystyle f(x)=\sum_i[a_i=x]\)，则
则 \(\displaystyle ans=\sum_{i=1}^n\sum_{j=1}^n\left\lfloor\dfrac{m}{a_ia_j}\right\rfloor = \sum_{a_ia_jk\le m}1 = \sum_{xyk\le m} f(x)f(y)\)。

其中 \(xyk\le m\) 意味着 \(x,y,k\) 中至少有一者 \(\le\sqrt[3]m\)。

记 \(B=\sqrt[3]m\)，作一下容斥：

\(\begin{aligned} ans &= 2\sum_{x\le B}f(x)\sum_{yk\le m/x}f(y) + \sum_{k\le B}\sum_{xy\le m/k}f(x)f(y) \\ &- 2\sum_{x\le B}f(x)\sum_{k\le B}Sf(m/(xk))-\sum_{x\le B}\sum_{y\le B}f(x)f(y)(m/(xy)) \\ &+ Sf(B)\cdot Sf(B)\cdot B. \end{aligned}\)

对于第一行的求和式，用整除分块计算、或者再容斥一次都可以，复杂度为 \(\displaystyle T(n)=\sum_{x\le B}\sqrt{\dfrac mx}=O(\sqrt{mB})=O(m^{2/3}).\)

需要用到 \(Sf(1)\sim Sf(\sqrt m),Sf(m/\sqrt m)\sim Sf(m/1)\) 的值（称为 \(f\) 的块筛），用类似桶排序的方法能在 \(O(n+\sqrt m)\) 的复杂度内求出。

因此总复杂度为 \(O(n+m^{2/3})\)。

#include<bits/stdc++.h>
#define rep(i, l, r) for (int i = l; i <= r; ++i)
typedef long long i64;

const int M = 1e5 + 8, P = 998244353;
int n, sq, B, x;
double inv[M];
i64 m, sf[M], Sf[M], sI[M], SI[M];

i64 dv(i64 x, int y) { return x * inv[y]; }

void init() {
    scanf("%d%lld", &n, &m);
    sq = sqrtl(m);
    B = cbrtl(m);
    rep(i, 1, n) {
        scanf("%d", &x);
        if (x <= sq) ++sf[x];
        else ++Sf[m / x];
    }
    rep(i, 1, sq) sf[i] += sf[i - 1];
    Sf[sq] += sf[sq];
    for (int i = sq; i; --i) Sf[i - 1] += Sf[i];
    rep(i, 1, sq) inv[i] = (1 + 1e-14) / i;
    rep(i, 1, sq) sI[i] = i;
    rep(i, 1, sq) SI[i] = dv(m, i);
}

i64 calc(i64* sf, int lim, i64* sg, i64* Sg, i64* sh, i64* Sh) {
    i64 ans = 0;
    rep(i, 1, B) {
        i64 m2 = dv(m, i);
        int sq2 = sqrtl(m2);
        int tmp = dv(sq, i);
        i64 sum = 0;
        if (!lim) {
            rep(j, 1, tmp) {
                sum += (sg[j] - sg[j - 1]) * Sh[i * j];
                sum += (sh[j] - sh[j - 1]) * Sg[i * j];
            }
            rep(j, tmp + 1, sq2) {
                sum += (sg[j] - sg[j - 1]) * sh[dv(m2, j)];
                sum += (sh[j] - sh[j - 1]) * sg[dv(m2, j)];
            }
            sum = (sum - sg[sq2] * sh[sq2]) % P;
        }
        else {
            tmp = std::min(tmp, lim);
            rep(j, 1, tmp) sum += (sg[j] - sg[j - 1]) * Sh[i * j];
            rep(j, tmp + 1, lim) sum += (sg[j] - sg[j - 1]) * sh[dv(m2, j)];
            sum %= P;
        }
        ans = (ans + (sf[i] - sf[i - 1]) * sum) % P;
    }
    return ans;
}

int main() {
    init();
    i64 ans = 0;
    ans += calc(sf, 0, sI, SI, sf, Sf) * 2;
    ans += calc(sI, 0, sf, Sf, sf, Sf);
    ans -= calc(sf, B, sI, SI, sf, Sf) * 2;
    ans -= calc(sf, B, sf, Sf, sI, SI);
    ans += B * sf[B] * sf[B];
    printf("%lld\n", (ans % P + P) % P);
    return 0;
}