自动搬运

查看原文

来自洛谷，原作者为

	XeCtera reverie

搬运于2025-08-24 22:49:04，当前版本为作者最后更新于2023-08-12 20:35:24，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

---

题目传送门：P9511 『STA - R3』大豆

0. 前言

给出一种时间复杂度 $O\left(\dfrac{km^{3/4}}{\log m}\right)$，空间复杂度 $O(\sqrt m)$ 的做法。

常数较小且不难实现。

1. 记号约定

记 $\displaystyle {\rm S}f(n)=\sum_{i=1}^nf(i)$。

记 $x/y=\lfloor\frac xy\rfloor$，且运算优先级低于乘法，即 $x/yz$ 表示 $x/(yz)$。

2. 题解 (Part 1)

不难发现 $\{a\}$ 的大豆化序列 $\{b\}$ 满足

$$b_m=a_m-\sum_{d=2}^mb_{m/d}$$

注意到这与杜教筛使用的式子具有相似性。

具体地，设 $I(n)=1$ 为常函数，$g*I=f$，则

$${\rm S}g(m)={\rm S}f(m)-\sum_{d=2}^m{\rm S}g(m/d)$$

当此式成立时，显然也能反推出 $g*I=f$。

本题中只需令 $f,g$ 分别为 $\{a\},\{b\}$ 的差分，那么有 $g=f*\mu$。

所求即为 ${\rm S}(f*\mu^k)(m)$ 的值。（$\mu^k$ 表示 $k$ 个 $\mu$ 卷积的结果）

使用足够厉害的数论科技即可解决。

下面开始介绍科技。

3. Dirichlet 双曲线法

这里可以把它当成整除分块的替代品。

相比于整除分块，它在常数等方面有一定优势。

考虑求如下式子的值：

$${\rm S}(f*g)(n)=\sum_{i=1}^nf(i){\rm S}g(n/i)=\sum_{ij\le n}f(i)g(j) $$

（假设 $f,g,{\rm S}f,{\rm S}g$ 的点值都可以 $O(1)$ 得到）

设阈值 $x,y$ 满足 $xy=n$。

我们对 $i\le x$ 的部分和 $j\le y$ 的部分分别求和，再减去公共部分。

$$\sum_{i\le x}f(i){\rm S}g(n/i)+\sum_{j\le y}g(j){\rm S}f(n/j)-{\rm S}f(x){\rm S}g(y) $$

可以 $O(x+y)$ 完成计算。

通常取 $x=y=\sqrt n$ 得最优复杂度 $O(\sqrt n)$。

4. 块筛卷积

首先我们有一条重要性质：$x/y/z=x/yz$。

根据这条性质，再结合整除分块的一些分析可知：

在求 ${\rm S}f(n)$ 时，我们往往只需用到 ${\rm S}f(1)\sim {\rm S}f(\lfloor\sqrt n\rfloor),{\rm S}f(n/\lfloor\sqrt n\rfloor)\sim{\rm S}f(n/1)$ 的值。

称之为 $f$（在 $n$ 处）的块筛。

块筛卷积：假设有一个狄利克雷卷积式 $f*g=h$，其中 $f,g$ 的块筛已知，试求 $h$ 的块筛。

可以直接用 Dirichlet 双曲线法，对块筛的每个位置暴力求解。

$${\rm S}h(n)=\sum_{i=1}^{\sqrt n}f(i){\rm S}g(n/i)+\sum_{i=1}^{\sqrt n}g(i){\rm S}f(n/i)-{\rm S}f(\sqrt n){\rm S}g(\sqrt n) $$

（开方运算默认下取整，下同）

复杂度是 $\displaystyle O\left(\sum_{i=1}^{\sqrt n}\sqrt i+\sum_{i=1}^{\sqrt n}\sqrt{n/i}\right)=O(n^{3/4})$。

可以进一步优化，但这里暂且不提。

5. 杜教筛

假设有一个狄利克雷卷积式 $f*g=h$，其中 $g,h$ 的块筛已知，试求 $f$ 的块筛。

事实上前面的做法仍然适用，只需要对式子做个移项：

$${\rm S}f(n)={\rm S}h(n)-\sum_{i=1}^{\sqrt n}f(i){\rm S}g(n/i)-\sum_{i=2}^{\sqrt n}g(i){\rm S}f(n/i)+{\rm S}f(\sqrt n){\rm S}g(\sqrt n) $$

递推即可，复杂度同样为 $O(n^{3/4})$。（同样可以优化，但先不提）

这种方法通常被称为杜教筛。

若使用普通数组存储块筛，并像上面这样递推求解，这种实现或许可以称为 dp 式杜教筛。

目前广为流传的则是递归 + map / unordered_map 记忆化存储的写法。我暂且称之为记搜式杜教筛。

事实上 dp 式实现的常数要小很多，并且复杂度也更加直观。

不明白为什么大家都用记搜式。难过。

6. 整除运算的优化

底层卡常技巧。

预处理 $1\sim\sqrt m$ 的倒数，将整除改为浮点乘。

若实现得当，可以保证代码中所有整除运算的除数都在 $1\sim\sqrt m$ 以内。将会有良好的优化效果。

注意事项：

用 (1+1e-15)/x 这样的方式来计算 x 的倒数，而非 1./x 。后者的精度损失可能导致出错。

7. 题解 (Part 2)

回到本题，我们要求 ${\rm S}(f*\mu^k)(m)$ 的值，其中 $f$ 的块筛已知。

朴素做法是对 $f$ 做 $k$ 次杜教筛。

事实上，dp 式杜教筛 + 整除优化的实现已经能够以 $O(km^{3/4})$ 的复杂度通过。

并且应该可以吊打 std。

但这还不够令人满意。下面给出一种更优秀的做法。

注意到 $\mu$ 是积性函数，比 $f$ 具有更好的性质。

考虑先求 $\mu^k$ 的块筛，然后容易 $O(\sqrt m)$ 得到 ${\rm S}(f*\mu^k)(m)$。

对于积性函数 $f$，记 $f_j(n)=f(n)[{\rm minp}(n)>p_j]$。

其中 ${\rm minp}(n)$ 是 $n$ 的最小质因子（钦定 ${\rm minp}(1)=+\infty$），$p_j$ 是第 $j$ 个质数。

取阈值 $K$ 满足 $p_{K+1}$ 略大于 $\sqrt[4]m$，那么 $f_K$ 有优秀的性质：

$f_K(1)\sim f_K(\sqrt m)$ 中，只有 $1$ 和大于 $\sqrt[4]m$ 的质数处的值非零。

回忆 Dirichlet 双曲线法的式子：

$${\rm S}h(n)=\sum_{i=1}^{\sqrt n}f(i){\rm S}g(n/i)+\sum_{i=1}^{\sqrt n}g(i){\rm S}f(n/i)-{\rm S}f(\sqrt n){\rm S}g(\sqrt n) $$

发现复杂度可优化至 $O(\pi(\sqrt n))=O\left(\dfrac{\sqrt n}{\ln n}\right)$。

我们用 $\mu_K*I_K=\varepsilon$ 杜教筛求出 $\mu_K$ 的块筛，再做块筛卷积得到 $\mu^k_K$ 的块筛。

这部分的复杂度为 $O\left(\dfrac{km^{3/4}}{\ln m}\right)$。

最后的问题是，如何得到 $I_K$ 的块筛，又如何将 $\mu_K^k$ 的块筛转化为 $\mu^k$ 的块筛？

事实上，采取 min_25 筛第一部分的 dp 即可实现。

给出 $I$ 的转移方程作为示例：

$$\begin{aligned} {\rm S}I_0(n)&={\rm S}I(n)\\&=n\\ {\rm S}I_j(n)&={\rm S}I_{j-1}(n)-I(p_j){\rm S}I_{j-1}(n/p_j)\\&={\rm S}I_{j-1}(n)-{\rm S}I_{j-1}(n/p_j)\quad(j>0) \end{aligned}$$

$\mu^k$ 同理，读者可以尝试自行推出。

这部分的复杂度为 $O(k\sqrt m\cdot\pi(\sqrt[4]m))=O\left(\dfrac{km^{3/4}}{\ln m}\right)$。

因此时间复杂度为 $O\left(\dfrac{km^{3/4}}{\ln m}\right)$。空间复杂度为 $O(\sqrt m)$。

8. 代码实现

变量名有些不太一致，注释以上文为准。

#include<bits/stdc++.h>
#define rep(i,l,r) for (int i=l; i<=r; ++i)
#define rrep(i,r,l) for (int i=r; i>=l; --i)
#define cpy(a,b) memcpy(a,b,sizeof a)
using namespace std;
typedef long long i64;

const int N=1e5+8,M=1e4+8,P=23068673;
const double I=1+1e-14;
i64 n; int T,K,sq,cnt,k,a[M],pri[M];
bitset<N/2>v; double inv[N];
int s1[N],s2[N],s3[N]; i64 S1[N],S2[N],S3[N];
struct Node { int f,g; }pos[M];

void init(int n) {
    int sq=sqrt(n),n2=(n>>1);
    for (int i=3; i<=sq; i+=2) if (!v[i>>1]) // 埃氏筛
        for (int j=(i*i>>1); j<=n2; j+=i) v[j]=1;
    pri[cnt=1]=2;
    rep(i,1,n2) if (!v[i]) pri[++cnt]=(i<<1|1);
    pri[cnt+1]=n+1;
    rep(i,1,n) inv[i]=I/i;
}
inline i64 dv(i64 x,int y) { return x*inv[y]; } // 整除优化

void calc1() {
    k=upper_bound(pri+1,pri+cnt+1,sqrt(sq))-pri;
    assert(pri[k+1]*pri[k+1]>sq);
    rep(i,1,sq) s1[i]=i;
    rep(i,1,sq) S1[i]=dv(n,i);
    rep(j,1,k) {
        int p=pri[j],t=dv(sq,p); i64 tn=dv(n,p);
        rep(i,1,t) S1[i]-=S1[i*p];
        rep(i,t+1,sq) S1[i]-=s1[dv(tn,i)];
        for (int i2=t,i=sq; i2; --i2)
            for (int lim=i2*p; i>=lim; --i) s1[i]-=s1[i2];
    } // 得到 I_K 的块筛
    rep(i,1,sq) s2[i]=2-s1[i];
    rrep(i,sq,1) {
        i64 n2=dv(n,i),s=0;
        int B=sqrt(n2),t,j=k+1;
        for ( ; (t=i*pri[j])<=sq; ++j)
            s+=S2[t]-S1[t];
        for ( ; pri[j]<=B; ++j)
            t=dv(n2,pri[j]),
            s+=s2[t]-s1[t];
        S2[i]=1-S1[i]-s+(i64)s1[B]*s2[B];
    } // 得到 μ_K 的块筛
    cpy(s1,s2),cpy(S1,S2);
}

void mul(int *sf,i64 *Sf,int *sg,i64 *Sg,int *sh,i64 *Sh) { // f*g=h 块筛卷积
    rep(j,k+1,cnt) {
        int p=pri[j];
        pos[j]={sf[p]-sf[p-1],sg[p]-sg[p-1]};
    }
    rep(i,1,sq) {
        i64 n2=dv(n,i),s=0;
        int B=sqrt(n2),t,j=k+1;
        for ( ; (t=i*pri[j])<=sq; ++j)
            s+=pos[j].f*Sg[t]+pos[j].g*Sf[t];
        for ( ; pri[j]<=B; ++j)
            t=dv(n2,pri[j]),
            s+=(i64)pos[j].f*sg[t]+(i64)pos[j].g*sf[t];
        Sh[i]=Sf[i]+Sg[i]+s-(i64)sf[B]*sg[B];
    }
    rrep(i,sq,1) sh[i]=sf[i]+sg[i]-1;
}

void solve() {
    rep(j,1,k) {
        int p=pri[j],t=dv(sq,p); i64 tn=dv(n,p);
        rep(tmp,1,K) {
            rep(i,1,t) S1[i]-=S1[i*p];
            rep(i,t+1,sq) S1[i]-=s1[dv(tn,i)];
            for (int i2=t,i=sq; i2; --i2)
                for (int lim=i2*p; i>=lim; --i) s1[i]-=s1[i2];
        }
    } // 得到 μ^k 的块筛
    rep(i,1,sq) s2[i]=a[(i-1)%T+1];
    rep(i,1,sq) S2[i]=a[(dv(n,i)-1)%T+1]; // 得到 f 的块筛
    i64 ans=0;
    rep(i,1,sq) ans=(ans+i64(s1[i]-s1[i-1])*S2[i])%P;
    rep(i,1,sq) ans=(ans+i64(s2[i]-s2[i-1])*S1[i])%P;
    ans=(ans-(i64)s1[sq]*s2[sq])%P;
    printf("%lld\n",ans<0?ans+P:ans);
}

int main() {
    scanf("%d%lld%d",&T,&n,&K);
    rep(i,1,T) scanf("%d",a+i),a[i]%=P;
    init(sq=sqrt(n));
    calc1();
    if (K>1) mul(s1,S1,s1,S1,s2,S2);
    if (K>2) mul(s1,S1,s2,S2,s3,S3);
    if (K==2) cpy(s1,s2),cpy(S1,S2);
    if (K==3) cpy(s1,s3),cpy(S1,S3);
    solve();
    return 0;
}