Min-25 筛学习笔记

Min-25 筛学习笔记

(text{By DaiRuiChen007})

一、简要介绍

Min-25 筛，是一种能在亚线性时间内求出特定的一类积性函数 (f(i)) 的前缀和的算法。

具体来说，Min-25 筛可以在 (mathcal O(sqrt n)) 的空间复杂度与 (mathcal O(dfrac{n^{3/4}}{log n})) 的时间复杂度内求出 (sum _{i=1}^n f(i)) 的值。

一般来说，积性函数 (f(i)) 要满足如下几个条件：

• (f(i)) 在质数处可以被表示成一个较为简单的多项式。
• (f(i)) 在质数的若干次幂（(p^c) 处）的点值可以 (mathcal O(1)) 计算。

一般来说，在 (n) 较小（(le 10^{10})）时，(mathcal O(dfrac{n^{3/4}}{log n})) 的复杂度是优于杜教筛的 (mathcal O(n^{2/3})) 的。

二、算法过程

以下内容均用 (p_i) 表示第 (i) 小的质数。

考虑拆分出每个 (n) 的最小质因子然后暴力计算，但是这样的质因子是 (mathcal O(n)) 级别的，考虑分成合数和质数两种情况讨论，其中合数的最小质因子就会被降到 (mathcal O(sqrt n)) 级别。

由于要根据最小质因子讨论，因此记 (S(n,k)) 表示 (1sim n) 中所有最小质因子 (>p_k) 的 (i) 的 (f(i)) 之和。

先不管质数的情况，设 (sumlimits_{i=1}^{p_ile n} f(p_i)) 的值为 (Q(n))，那么枚举 (>p_k) 的最小质因子以及其指数可以得到：

[S(n,k)=Q(n)-sum_{i=1}^{ile k} f(p_i)+sum_{i=k+1}^{p_i^2le n}sum_{c=1}^{p_i^cle n} f(p_i^c)times (S(lfloor n/p_i^crfloor,i)+[c>1]) ]

其中 ([c>1]) 是求 (f(p_i^c)) 的贡献，容易发现答案就是 (S(n,0)+f(1))。

假设 (Q(n)) 已知，预处理出 (f(p_i)) 的前缀和（容易递归过程中发现 (p_k^2le n)，因此只要处理 (le sqrt n) 的质数），这个式子直接暴力递归计算的复杂度就是 (mathcal O(dfrac{n^{3/4}}{log n}))。

接下来考虑计算 (Q(n))，注意到在递归过程中的 (n) 始终可以表示成某个 (lfloor n/trfloor) 的值，因此我们只要求出 (Q(n)) 在这 (2sqrt n) 个位置的点值即可。

根据要求，(f(p_i)) 一定是一个简单的多项式，我们分次数处理，分别求出 $sum p_i^0,sum p_i^1,sum p_i^2,dots $ 再整合。以 (sum p_i^1) 为例，构造完全积性函数 (f'(i)=i)，容易发现这个函数在质数处与 (f(i)) 相等，因此我们只要求出 (f'(i)) 在质数处的和。

设 (G(n,k)) 表示 (1sim n) 中，满足 (i) 的最小质因子 (>p_k) 或 (i) 是质数的 (f'(i)) 的和。

考虑 (G(n,k-1)-G(n,k)) 的值，当 (p_k^2>n) 时，显然没有最小质因子 (=p_k) 的合数，差为 (0)。

否则，计算的数一定有最小质因子 (p_k)，提出 (f'(p_k)) 得到：

[G(n,k-1)-G(n,k)=f'(p_k)timesleft(G(lfloor n/p_krfloor,k-1)-sum_{i=1}^{k-1}f'(p_i) right) ]

减去的项是去掉有 (<p_k) 的质因子的项，因此我们得到 (G(n,k)) 的递推式：

[G(n,k)= begin{cases} G(n,k-1)&p_k^2>n\ G(n,k-1)-f'(p_k)timesleft(G(lfloor n/p_krfloor,k-1)-sum_{i=1}^{k-1}f'(p_i) right) &p_k^2le n end{cases} ]

同样在递归的过程中也只会访问到所有的 (lfloor n/trfloor) 下标处，且最后 (f'(p_i)) 的前缀和也只要处理到 (le sqrt n) 的所有质数。

边界条件是 (G(n,0)=n(n+1)/2-1)，(Q(n)=G(n,infty))，从小到大枚举 (k)，从大到小枚举 (n) 转移即可，复杂度也可以证明为 (mathcal O(dfrac{n^{3/4}}{log n}))。

拓展 - 非递归版本的实现

对于 (f(i)) 的前缀和 (S(n))，我们已经介绍了在 (mathcal O(dfrac{n^{3/4}}{log n})) 的时间复杂度以及 (mathcal O(sqrt n)) 的空间复杂度内计算 (S(n)) 的算法。

但假如我们想和杜教筛一样，在复杂度不变的情况下求出 (S(n)) 在所有 (S(lfloor n/trfloor )) 处的值，应该怎么处理呢？

最简单的想法是在递归求 (S(n,k)) 的时候加上记忆化，这样确实能保证时间复杂度不退化，但是空间复杂度会变成 (mathcal O(dfrac{n^{3/4}}{log n}))，如果想要空间复杂度是 (mathcal O(sqrt n)) 的做法，需要转变一下 (S) 的定义再处理。

注意到在处理 (S(n,k)) 时最棘手的是形如 (S(n,k)gets S(n',i)) 的转移时 (k) 与 (i) 不连续，因此导致没法简单递推求值。

我们认为 (G) 的递推式很优秀也是因为转移式永远只有 (i=k-1) 的状态向 (k) 转移，因此考虑让 (S(n,k)) 的定义也向 (G(n,k)) 的定义靠拢，不妨设 (S(n,k)) 表示 (1sim n) 中，满足 (i) 的最小质因子 (>p_k) 或 (i) 是质数的 (f(i)) 的和。

依然考虑 (S(n,k-1)-S(n,k))，在 (p_k^2> n) 时依然为 (0)，在 (p_k^2le n) 时枚举 (p_k) 的次数得到：

[S(n,k-1)-S(n,k)=sum_{c=1}^{p_k^cle n} f(p_k^c)timesleft(S(lfloor n/p_k^crfloor,k)+[c>1]-sum_{i=1}^{p_ile min(p_k,lfloor n/p_k^crfloor)}f(p_i)right) ]

但是注意到最后一项求和式并不是很好看，考虑化简，注意到当 (p_k> lfloor n/p_k^crfloor) 时，(p_k^{c+1}> n)，那么 (S(lfloor n/p_k^crfloor,k)) 只会统计 (le lfloor n/p_k^crfloor) 的质数，与求和式正好抵消，因此可以把上式写成：

[S(n,k-1)-S(n,k)=sum_{c=1}^{p_k^{c+1}le n}f(p_k^{c+1})+ f(p_k^c)timesleft(S(lfloor n/p_k^crfloor,k)-sum_{i=1}^{k}f(p_i)right) ]

因此得到 (S(n,k)) 的递推式：

[S(n,k-1)= begin{cases} S(n,k)&p_k^2>n\[2ex] S(n,k)+sumlimits_{c=1}^{p_k^{c+1}le n}f(p_k^{c+1})+ f(p_k^c)timesleft(S(lfloor n/p_k^crfloor,k)-sum_{i=1}^{k}f(p_i)right)&p_k^2le n end{cases} ]

从大到小枚举 (p_k) 再从大到小枚举 (n) 转移即可。

其中边界条件为 (S(n,infty)=Q(n))，最终的答案是 (S(lfloor n/trfloor,0))，时间复杂度 (mathcal O(dfrac{n^{3/4}}{log n}))，空间复杂度 (mathcal O(sqrt n))。

有了这个 Min-25 筛的实现后，我们可以解决形如 (sum_{i=1}^n f(i)times g(lfloor n/irfloor)) 的和式计算问题，其中 (f(i)) 是积性函数，(g(i)) 是任意一个关于 (i) 的函数。

那么我们根据整除分块枚举 $lfloor n/irfloor $，则 (g) 函数的值固定，要求的就是某一段特定区间内 (f(i)) 的和，容易发现这些区间的端点都是某些 (lfloor n/trfloor)，用上面所介绍的非递归型的 Min-25 筛处理即可。

三、代码实现

以线性筛模板题为例（Link）下面的代码是递归的实现：

#include<bits/stdc++.h> #define int long long using namespace std; const int MAXN=2e5+1,MOD=1e9+7; bool isco[MAXN]; int p[MAXN],tot; int val[MAXN]; //需要处理的 n (降序排列) int idx1[MAXN],idx2[MAXN]; //用于储存每个可能的 n/t 对应的下标 int g1[MAXN],g2[MAXN]; //p,p^2 在所有 n/t 处的前缀和 int fp1[MAXN],fp2[MAXN]; //p,p^2 在所有 <=sqrt(n) 的质数处的前缀和 int n,m,B; inline int idx(int v) {     //v 对应存储在数组里的下标     return (v<=B)?idx1[v]:idx2[n/v]; } inline int S(int n,int k) {     if(p[k]>n) return 0;     int ans=((g2[idx(n)]+MOD-g1[idx(n)])%MOD+MOD-(fp2[k]+MOD-fp1[k])%MOD)%MOD;     for(int i=k+1;i<=tot&&p[i]*p[i]<=n;++i) {         for(int c=1,v=p[i];v<=n;++c,v*=p[i]) {             //直接递归求值             ans=(ans+(v-1)%MOD*(v%MOD)%MOD*((c>1)+S(n/v,i))%MOD)%MOD;         }     }     return ans; } signed main() {     scanf("%lld",&n),B=sqrt(n);     for(int i=2;i<=B;++i) {         if(!isco[i]) p[++tot]=i;         for(int j=1;j<=tot&&i*p[j]<=B;++j) {             isco[i*p[j]]=true;             if(i%p[j]==0) break;         }         //线性筛预处理质数     }     for(int i=1;i<=tot;++i) {         //在 <=sqrt(n) 的质数处的前缀和         fp1[i]=(fp1[i-1]+p[i])%MOD;         fp2[i]=(fp2[i-1]+p[i]*p[i]%MOD)%MOD;     }     for(int l=1,r;l<=n;l=r+1) {         //整除分块预处理出点值         r=n/(n/l),val[++m]=n/l;         if(val[m]<=B) idx1[val[m]]=m;         else idx2[n/val[m]]=m;     }     for(int i=1;i<=m;++i) {         int k=val[i]%MOD;         g1[i]=(MOD+1)/2*k%MOD*(k+1)%MOD;         g2[i]=(MOD+1)/6*k%MOD*(2*k+1)%MOD*(k+1)%MOD;         g1[i]=g1[i]?g1[i]-1:MOD-1,g2[i]=g2[i]?g2[i]-1:MOD-1;         //G(n,0) 的值     }     for(int k=1;k<=tot;++k) {         for(int i=1;i<=m&&p[k]*p[k]<=val[i];++i) {             //暴力枚举, 按递推式计算             g1[i]=(g1[i]+MOD-p[k]*(g1[idx(val[i]/p[k])]+MOD-fp1[k-1])%MOD)%MOD;             g2[i]=(g2[i]+MOD-p[k]*p[k]%MOD*(g2[idx(val[i]/p[k])]+MOD-fp2[k-1])%MOD)%MOD;         }     }     printf("%lldn",(S(n,0)+1)%MOD);     return 0; } 

下面的代码是非递归的实现：

#include<bits/stdc++.h> #define int long long using namespace std; const int MAXN=2e5+1,MOD=1e9+7; bool isco[MAXN]; int p[MAXN],tot; int val[MAXN];  int idx1[MAXN],idx2[MAXN]; int g1[MAXN],g2[MAXN]; int fp1[MAXN],fp2[MAXN]; int S[MAXN]; int n,m,B; inline int idx(int v) {     return (v<=B)?idx1[v]:idx2[n/v]; } signed main() {     scanf("%lld",&n),B=sqrt(n);     for(int i=2;i<=B;++i) {         if(!isco[i]) p[++tot]=i;         for(int j=1;j<=tot&&i*p[j]<=B;++j) {             isco[i*p[j]]=true;             if(i%p[j]==0) break;         }     }     for(int i=1;i<=tot;++i) {         fp1[i]=(fp1[i-1]+p[i])%MOD;         fp2[i]=(fp2[i-1]+p[i]*p[i]%MOD)%MOD;     }     for(int l=1,r;l<=n;l=r+1) {         r=n/(n/l),val[++m]=n/l;         if(val[m]<=B) idx1[val[m]]=m;         else idx2[n/val[m]]=m;     }     for(int i=1;i<=m;++i) {         int k=val[i]%MOD;         g1[i]=(MOD+1)/2*k%MOD*(k+1)%MOD;         g2[i]=(MOD+1)/6*k%MOD*(2*k+1)%MOD*(k+1)%MOD;         g1[i]=g1[i]?g1[i]-1:MOD-1,g2[i]=g2[i]?g2[i]-1:MOD-1;     }     for(int k=1;k<=tot;++k) {         for(int i=1;i<=m&&p[k]*p[k]<=val[i];++i) {             g1[i]=(g1[i]+MOD-p[k]*(g1[idx(val[i]/p[k])]+MOD-fp1[k-1])%MOD)%MOD;             g2[i]=(g2[i]+MOD-p[k]*p[k]%MOD*(g2[idx(val[i]/p[k])]+MOD-fp2[k-1])%MOD)%MOD;         }     }     for(int i=1;i<=m;++i) S[i]=(g2[i]+MOD-g1[i])%MOD;     for(int k=tot;k>=1;--k) {         for(int i=1;i<=m&&p[k]*p[k]<=val[i];++i) {             for(int c=1,v=p[k];v*p[k]<=val[i];) { //根据递推式计算                 S[i]=(S[i]+(v-1)%MOD*(v%MOD)%MOD*(S[idx(val[i]/v)]+MOD-(fp2[k]-fp1[k]))%MOD)%MOD;                 ++c,v*=p[k],S[i]=(S[i]+(v-1)%MOD*(v%MOD)%MOD)%MOD;             }         }     }     printf("%lldn",(S[idx(n)]+1)%MOD);     return 0; } 

实际运行上，非递归版本的实现会慢于递归版的实现。

四、复杂度分析

空间复杂度 (mathcal O(sqrt n))，可以参照参考代码的实现。

时间复杂度可以看做合法的 (S(n,k) / G(n,k)) 状态数，对于每个 (n=lfloor n/trfloor)，合法的 (k) 共有 (pi(sqrt{lfloor n/trfloor})) 个。

因此复杂度可以估计为：

[begin{aligned} mathcal T(n) &=sum_{i^2le n} mathcal O(pi(sqrt i))+mathcal O(pi(sqrt{n/i}))\[2.5ex] &=sum_{i^2le n} mathcal Oleft(dfrac{sqrt i}{ln sqrt i}right)+mathcal Oleft(dfrac{sqrt{n/i}}{ln sqrt{n/i}}right)\[2.5ex] &=mathcal Oleft(dfrac{1}{log n}int_1^{sqrt n} sqrt{dfrac nx} mathrm dxright)\[2.5ex] &=mathcal Oleft(dfrac{sqrt n}{log n}timesleft. 2sqrt x right|^{sqrt n}_1 right)\[2.5ex] &=mathcal O(dfrac{n^{3/4}}{log n}) end{aligned} ]

五、例题

I. [SPOJ-34096] DIVCNTK

Problem Link

题目大意

给定 (n)，求 (sum_{i=1}^n sigma_0(i^k))（即 (i^k) 的因子个数）。

数据范围：(nle 10^{10})。

思路分析

容易证明 (sigma_0(i)) 是积性函数，同理 (sigma_0(i^k)) 同样也是积性函数，考虑求 (sigma_0(p^k))，显然 (sigma_0(p^k)=k+1)，我们只需要求所有 (1sim lfloor n/trfloor) 之间的质数个数，构造 (f'(p)=1)，直接用 Min-25 筛处理即可。

时间复杂度 (mathcal O(dfrac{n^{3/4}}{log n}))

代码呈现

#include<bits/stdc++.h> #define int unsigned long long using namespace std; const int MAXN=2e5+1; bool isco[MAXN]; int p[MAXN],tot; int idx1[MAXN],idx2[MAXN],val[MAXN]; int g[MAXN]; inline void solve() {     int n,K,m=0,B;     scanf("%llu%llu",&n,&K),B=sqrt(n);     for(int l=1,r;l<=n;l=r+1) {         r=n/(n/l),val[++m]=n/l;         if(val[m]<=B) idx1[val[m]]=m;         else idx2[n/val[m]]=m;         g[m]=val[m]-1;     }     auto idx=[&](int v) { return v<=B?idx1[v]:idx2[n/v]; };     for(int k=1;k<=tot;++k) {         for(int i=1;i<=m&&p[k]*p[k]<=val[i];++i) {             g[i]-=g[idx(val[i]/p[k])]-(k-1);         }     }     auto S=[&](auto self,int n,int k) -> int {         if(n<=p[k]) return 0;         int ans=(g[idx(n)]-k)*(K+1);         for(int i=k+1;i<=tot&&p[i]*p[i]<=n;++i) {             for(int c=1,v=p[i];v<=n;++c,v*=p[i]) {                 ans+=(c*K+1)*(self(self,n/v,i)+(c>1));             }         }         return ans;     };     printf("%llun",S(S,n,0)+1); } signed main() {     for(int i=2;i<MAXN;++i) {         if(!isco[i]) p[++tot]=i;         for(int j=1;j<=tot&&i*p[j]<MAXN;++j) {             isco[i*p[j]]=true;             if(i%p[j]==0) break;         }     }     int T;     scanf("%llu",&T);     while(T--) solve();     return 0; } 

---

II. [洛谷-4213] SUM

Problem Link

题目大意

求 (mu(i),varphi(i)) 的前缀和。

数据范围：(nle 2^{31})。

思路分析

注意到 (mu(p^c),varphi(p^c)) 的计算是容易的，在质数处的点值满足：(mu(p)=-1,varphi(p)=p-1)，因此预处理 (G(n,k)) 时分别处理 (sum p_i^0,sum p_i^1) 的值即可。

时间复杂度 (mathcal O(dfrac{n^{3/4}}{log n}))。

代码呈现

#include<bits/stdc++.h> #define int long long using namespace std; const int MAXN=2e5+1; bool isco[MAXN]; int p[MAXN],tot,val[MAXN],fp1[MAXN]; int idx1[MAXN],idx2[MAXN],g0[MAXN],g1[MAXN]; inline void solve() {     int n,m=0;     scanf("%lld",&n);     int B=sqrt(n); tot=0;     for(int i=2;i<=B;++i) {         if(!isco[i]) p[++tot]=i,fp1[tot]=i+fp1[tot-1];         for(int j=1;j<=tot&&p[j]*i<=B;++j) {             isco[p[j]*i]=true;             if(i%p[j]==0) break;         }     }     for(int l=1,r;l<=n;l=r+1) {         r=n/(n/l),val[++m]=n/l;         if(val[m]<=B) idx1[val[m]]=m;         else idx2[n/val[m]]=m;         g0[m]=val[m]-1,g1[m]=val[m]*(val[m]+1)/2-1;     }     auto idx=[&](int v) { return (v<=B)?idx1[v]:idx2[n/v]; };     for(int k=1;k<=tot;++k) {         for(int i=1;i<=m&&p[k]*p[k]<=val[i];++i) {             g0[i]-=g0[idx(val[i]/p[k])]-(k-1);             g1[i]-=p[k]*(g1[idx(val[i]/p[k])]-fp1[k-1]);         }     }     auto mu=[&](auto self,int n,int k) -> int {         if(n<=p[k]) return 0;         int ans=k-g0[idx(n)];         for(int i=k+1;i<=tot&&p[i]*p[i]<=n;++i) {             ans-=self(self,n/p[i],i);         }         return ans;     };     auto phi=[&](auto self,int n,int k) -> int {         if(n<=p[k]) return 0;         int ans=g1[idx(n)]-g0[idx(n)]-(fp1[k]-k);         for(int i=k+1;i<=tot&&p[i]*p[i]<=n;++i) {             for(int c=1,v=p[i];v<=n;++c,v*=p[i]) {                 ans+=(v-v/p[i])*(self(self,n/v,i)+(c>1));             }         }         return ans;     };     printf("%lld %lldn",phi(phi,n,0)+1,mu(mu,n,0)+1); } signed main() {     int T;     scanf("%lld",&T);     while(T--) solve();     return 0; } 

---

* III. [HDU-6417] Rikka with APSP

Problem Link

题目大意

给定一张 (n) 个点的无向完全图，(i,j) 之间的边权定义为最小的正整数 (x) 使得 (ijx) 是完全平方数。

求 (sum_{1le ile jle n}mathrm{dist}(i,j)) 的值，其中 (mathrm{dist}(i,j)) 表示 ((i,j)) 之间的最短路长度。

数据范围：(nle 10^{10})。

思路分析

先把每个数 (k) 写成标准分解形式 (prod p_i^{c_i})，注意到我们只关心 (c_ibmod 2) 的值，因此可以把所有的这些 (c_ibmod 2) 写成一个二进制数 (B_k)。

那么 (w(u,v)) 就是 (B_uoplus B_v) 中所有 (1) 对应的 (p_i) 的乘积，而 (mathrm{dist}(u,v)) 就是 (B_uoplus B_v) 中所有 (1) 对应的 (p_i) 的乘积，先逐步去掉 (B_ucap(B_uoplus B_v)) 中的 (1)，再逐步加入 (B_vcap(B_uoplus B_v)) 中的 (1)，容易证明整个过程中的数值 (le max(u,v))。

因此考虑拆贡献，对于某个质数 (p)，记 (f(p)) 为 ([1,n]) 中含有奇数个 (p) 的数的数量，那么 (p) 对答案的贡献就是 (ptimes f(p)times(n-f(p)))。

容易发现 (f(p)=leftlfloordfrac{n}{p}rightrfloor-leftlfloordfrac{n}{p^2}rightrfloor+leftlfloordfrac{n}{p^3}rightrfloor-dots)，考虑根号分治，对于 (ple sqrt n) 的 (p)，暴力计算，时间复杂度为 (mathcal O(pi(sqrt n)log n)=mathcal O(dfrac{sqrt n}{log n}log n)=mathcal O(sqrt n))。

对于 (p>sqrt n) 的 (p)，(f(p)=leftlfloordfrac nprightrfloor)，用整除分块的技巧枚举 (leftlfloordfrac nprightrfloor=k)，那么 (f(p)times (n-f(p))) 容易计算，那么我们只要求某个特定区间 ([l,r]) 之中的质数之和，考虑整除分块的过程，注意到 (l-1,r) 都能写成某个 (leftlfloordfrac ntrightrfloor) 的形式。

因此我们只预处理在所有 (1sim leftlfloordfrac ntrightrfloor) 之间质数之和，在整除分块的过程中就可以 (mathcal O(1)) 回答对应的区间素数数量了。

容易发现这是一个标准的 Min-25 筛的第一部分，因此直接用 Min-25 筛处理即可。

这一部分的时间复杂度为 (mathcal O(dfrac{n^{3/4}}{log n}))。

但是我们还漏掉了一些贡献，对于 (utimes v) 是完全平方数但 (une v) 的数对 ((u,v)) 对答案是有 (1) 的贡献的，因此我们要求，有多少对 ((u,v)) 满足 (1le u<vle n)，且 (utimes v) 是完全平方数，考虑构造一个等比数列 ({v,sqrt{uv},u})，其公比 (<1) 且所有元素都是 ([1,n]) 之间正整数，我们只需要对这样的等比数列计数即可。

设公比最简分数为 (q/p)，枚举分母 (p)，此时分子 (q<p) 且 (gcd(p,q)=1)，因此合法的 (q) 恰好有 (varphi(p)) 种选择。

然后考虑 (v)，注意到 (u=dfrac {vq^2}{p^2})，因此 (p^2mid v)，这样的 (v) 共有 (leftlfloordfrac n{p^2}rightrfloor) 种，容易证明这样的每一对 ((p,q,v)) 都和我们要计数的 ((u,v)) 构成双射。

因此这部分的答案就是 (sum_{p=2}^{sqrt n} varphi(p)times leftlfloordfrac n{p^2}rightrfloor)，预处理出 (varphi(1)sim varphi(sqrt n)) 即可做到 (mathcal O(sqrt n))。

将答案分成三部分分别处理即可。

第三部分的另一种处理方法

对于这样乘积是完全平方数的数对计数，最朴素的思路就是枚举 (u,v) 在各自除掉最大的完全平方因子后得到的数 (k)，枚举这样的 (k)，(u,v) 的选择方案数就是 (g(lfloor n/krfloor)=binom{sqrt{lfloor n/krfloor}}2)。

考虑刻画 (k) 以写出求和式，显然 (k) 中无平方因子，则 (mu(k)ne 0)，因此满足条件的 (k) 可以用 (mu^2(k)=1) 表示，因此这部分的贡献就是：

[sum_{k=1}^n mu^2(k)binom{sqrt{lfloor n/krfloor}}2 ]

考虑整除分块枚举 (lfloor n/krfloor)，原问题变成统计 (mu^2(k)) 在所有 (lfloor n/trfloor) 处的前缀和，注意到 (mu^2(k)) 是积性函数，因此可以用非递归版本的 Min-25 筛求出。

由于 (mu^2(k)) 函数的特殊性，在递推 (S(n,k)) 的时候只需要处理 (c=1) 的情况，因此该算法的运行效率和前一个算法的运行效率差别并不大。

两种做法时间复杂度均为 (mathcal O(dfrac{n^{3/4}}{log n}))。

代码呈现

实现方式一：

#include<bits/stdc++.h> #define int long long #define LL __int128 using namespace std; const int MAXN=2e5+1,MOD=998244353; bool isco[MAXN]; int p[MAXN],tot,val[MAXN],fp[MAXN],phi[MAXN]; int idx1[MAXN],idx2[MAXN],g[MAXN]; inline void solve() {     int n,m=0;     scanf("%lld",&n);     int B=sqrt(n);     for(int l=1,r;l<=n;l=r+1) {         r=n/(n/l),val[++m]=n/l;         if(val[m]<=B) idx1[val[m]]=m;         else idx2[n/val[m]]=m;     }     auto idx=[&](int v) { return (v<=B)?idx1[v]:idx2[n/v]; };     for(int i=1;i<=m;++i) {         int k=val[i]%MOD;         g[i]=(k*(k+1)/2+MOD-1)%MOD;     }     for(int k=1;k<=tot;++k) {         for(int i=1;i<=m&&p[k]*p[k]<=val[i];++i) {             g[i]=(g[i]+MOD-p[k]*(g[idx(val[i]/p[k])]+MOD-fp[k-1])%MOD)%MOD;         }     }     int ans=0;     for(int i=1;i<=tot&&p[i]*p[i]<=n;++i) {         int s=0;         for(int v=p[i],c=1;v<=n;v*=p[i],c=-c) s+=c*(n/v);         ans=(ans+(n-s)%MOD*(s%MOD)%MOD*p[i]%MOD)%MOD;     }     for(int l=1,r;l<=n;l=r+1) {         r=n/(n/l); int k=n/l;         if((LL)l*l>n) {             ans=(ans+(n-k)%MOD*(k%MOD)%MOD*(g[idx(r)]+MOD-g[idx(l-1)])%MOD)%MOD;             }     }     for(int i=2;i*i<=n;++i) ans=(ans+(n/(i*i))%MOD*phi[i]%MOD)%MOD;     printf("%lldn",ans); } signed main() {     for(int i=2;i<MAXN;++i) {         if(!isco[i]) p[++tot]=i,phi[i]=i-1;         for(int j=1;j<=tot&&p[j]*i<MAXN;++j) {             isco[p[j]*i]=true,phi[p[j]*i]=phi[i]*(p[j]-1);             if(i%p[j]==0) { phi[p[j]*i]=phi[i]*p[j]; break; }         }     }     for(int i=1;i<=tot;++i) fp[i]=(fp[i-1]+p[i])%MOD;     int T;     scanf("%lld",&T);     while(T--) solve();     return 0; } 

实现方式二：

#include<bits/stdc++.h> #define int long long #define LL __int128 using namespace std; const int MAXN=2e5+1,MOD=998244353; bool isco[MAXN]; int p[MAXN],tot,val[MAXN],fp1[MAXN]; int idx1[MAXN],idx2[MAXN],g0[MAXN],g1[MAXN],S[MAXN]; inline void solve() {     int n,m=0;     scanf("%lld",&n);     int B=sqrt(n); tot=0;     for(int i=2;i<=B;++i) {         if(!isco[i]) p[++tot]=i;         for(int j=1;j<=tot&&p[j]*i<=B;++j) {             isco[p[j]*i]=true;             if(i%p[j]==0) break;         }     }     for(int i=1;i<=tot;++i) fp1[i]=(fp1[i-1]+p[i])%MOD;     for(int l=1,r;l<=n;l=r+1) {         r=n/(n/l),val[++m]=n/l;         if(val[m]<=B) idx1[val[m]]=m;         else idx2[n/val[m]]=m;         int k=val[m]%MOD;         g0[m]=k-1,g1[m]=(k*(k+1)/2+MOD-1)%MOD;     }     auto idx=[&](int v) { return (v<=B)?idx1[v]:idx2[n/v]; };     for(int k=1;k<=tot;++k) {         for(int i=1;i<=m&&p[k]*p[k]<=val[i];++i) {             g0[i]=(g0[i]+MOD-(g0[idx(val[i]/p[k])]-(k-1)))%MOD;             g1[i]=(g1[i]+MOD-p[k]*(g1[idx(val[i]/p[k])]+MOD-fp1[k-1])%MOD)%MOD;         }     }     for(int i=1;i<=m;++i) S[i]=g0[i];     for(int k=tot;k>=1;--k) {         for(int i=1;i<=m&&p[k]*p[k]<=val[i];++i) {             S[i]=(S[i]+MOD+S[idx(val[i]/p[k])]-k)%MOD;         }     }     for(int i=1;i<=m;++i) S[i]=(S[i]+1)%MOD;     int ans=0;     for(int i=1;i<=tot&&p[i]*p[i]<=n;++i) {         int s=0;         for(int v=p[i],c=1;v<=n;v*=p[i],c=-c) s+=c*(n/v);         ans=(ans+(n-s)%MOD*(s%MOD)%MOD*p[i]%MOD)%MOD;     }     for(int l=1,r;l<=n;l=r+1) {         r=n/(n/l); int k=n/l;         if((LL)l*l>n) {             ans=(ans+(n-k)%MOD*(k%MOD)%MOD*(g1[idx(r)]+MOD-g1[idx(l-1)])%MOD)%MOD;             }     }     for(int l=1,Q=B+5,r;l<=n;l=r+1) {         r=n/(n/l); int k=n/l;         while(Q*Q>k) --Q; //Q = sqrt(k)         ans=(ans+(Q*(Q-1)/2)%MOD*(S[idx(r)]+MOD-S[idx(l-1)])%MOD)%MOD;     }     printf("%lldn",ans); } signed main() {     int T;     scanf("%lld",&T);     while(T--) solve();     return 0; }