Min-25 筛学习笔记

$By DaiRuiChen007$

一、简要介绍

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

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

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

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

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

二、算法过程

以下内容均用 $p_{i}$ 表示第 $i$ 小的质数。

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

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

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

S (n, k) = Q (n) - \sum_{i = 1}^{i \leq k} f (p_{i}) + \sum_{i = k + 1}^{p_{i}^{2} \leq n} \sum_{c = 1}^{p_{i}^{c} \leq n} f (p_{i}^{c}) \times (S (⌊ n / p_{i}^{c} ⌋, i) + [c > 1])

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

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

接下来考虑计算 $Q (n)$ ，注意到在递归过程中的 $n$ 始终可以表示成某个 $⌊ n / t ⌋$ 的值，因此我们只要求出 $Q (n)$ 在这 $2 \sqrt{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)$ 表示 $1 \sim 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}) \times (G (⌊ n / p_{k} ⌋, k - 1) - \sum_{i = 1}^{k - 1} f^{'} (p_{i}))

减去的项是去掉有 $< 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}) \times (G (⌊ n / p_{k} ⌋, k - 1) - \sum_{i = 1}^{k - 1} f^{'} (p_{i})) & p_{k}^{2} \leq n \end{cases}

同样在递归的过程中也只会访问到所有的 $⌊ n / t ⌋$ 下标处，且最后 $f^{'} (p_{i})$ 的前缀和也只要处理到 $\leq \sqrt{n}$ 的所有质数。

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

拓展 – 非递归版本的实现

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

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

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

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

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

依然考虑 $S (n, k - 1) - S (n, k)$ ，在 $p_{k}^{2} > n$ 时依然为 $0$ ，在 $p_{k}^{2} \leq n$ 时枚举 $p_{k}$ 的次数得到：

S (n, k - 1) - S (n, k) = \sum_{c = 1}^{p_{k}^{c} \leq n} f (p_{k}^{c}) \times (S (⌊ n / p_{k}^{c} ⌋, k) + [c > 1] - \sum_{i = 1}^{p_{i} \leq min (p_{k}, ⌊ n / p_{k}^{c} ⌋)} f (p_{i}))

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

S (n, k - 1) - S (n, k) = \sum_{c = 1}^{p_{k}^{c + 1} \leq n} f (p_{k}^{c + 1}) + f (p_{k}^{c}) \times (S (⌊ n / p_{k}^{c} ⌋, k) - \sum_{i = 1}^{k} f (p_{i}))

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

S (n, k - 1) = {\begin{cases} S (n, k) & p_{k}^{2} > n \\ S (n, k) + \sum_{c = 1}^{p_{k}^{c + 1} \leq n} f (p_{k}^{c + 1}) + f (p_{k}^{c}) \times (S (⌊ n / p_{k}^{c} ⌋, k) - \sum_{i = 1}^{k} f (p_{i})) & p_{k}^{2} \leq n \end{cases}

从大到小枚举 $p_{k}$ 再从大到小枚举 $n$ 转移即可。

其中边界条件为 $S (n, \infty) = Q (n)$ ，最终的答案是 $S (⌊ n / t ⌋, 0)$ ，时间复杂度 $O (\frac{n^{3 / 4}}{\log n})$ ，空间复杂度 $O (\sqrt{n})$ 。

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

那么我们根据整除分块枚举 $⌊ n / i ⌋$ ，则 $g$ 函数的值固定，要求的就是某一段特定区间内 $f (i)$ 的和，容易发现这些区间的端点都是某些 $⌊ n / t ⌋$ ，用上面所介绍的非递归型的 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("%lld\n",(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("%lld\n",(S[idx(n)]+1)%MOD);
return 0;
}

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

四、复杂度分析

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

时间复杂度可以看做合法的 $S (n, k) / G (n, k)$ 状态数，对于每个 $n = ⌊ n / t ⌋$ ，合法的 $k$ 共有 $π (\sqrt{⌊ n / t ⌋})$ 个。

因此复杂度可以估计为：

\begin{aligned} T (n) & = \sum_{i^{2} \leq n} O (π (\sqrt{i})) + O (π (\sqrt{n / i})) \\ = \sum_{i^{2} \leq n} O (\frac{\sqrt{i}}{\ln \sqrt{i}}) + O (\frac{\sqrt{n / i}}{\ln \sqrt{n / i}}) \\ = O (\frac{1}{\log n} \int_{1}^{\sqrt{n}} \sqrt{\frac{n}{x}} d x) \\ = O (\frac{\sqrt{n}}{\log n} \times {2 \sqrt{x} |}_{1}^{\sqrt{n}}) \\ = O (\frac{n^{3 / 4}}{\log n}) \end{aligned}

五、例题

I. [SPOJ-34096] DIVCNTK

Problem Link

题目大意

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

数据范围： $n \leq 10^{10}$ 。

思路分析

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

时间复杂度 $O (\frac{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("%llu\n",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

题目大意

求 $μ (i), φ (i)$ 的前缀和。

数据范围： $n \leq 2^{31}$ 。

思路分析

注意到 $μ (p^{c}), φ (p^{c})$ 的计算是容易的，在质数处的点值满足： $μ (p) = - 1, φ (p) = p - 1$ ，因此预处理 $G (n, k)$ 时分别处理 $\sum p_{i}^{0}, \sum p_{i}^{1}$ 的值即可。

时间复杂度 $O (\frac{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 %lld\n",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$ 使得 $i j x$ 是完全平方数。

求 $\sum_{1 \leq i \leq j \leq n} dist (i, j)$ 的值，其中 $dist (i, j)$ 表示 $(i, j)$ 之间的最短路长度。

数据范围： $n \leq 10^{10}$ 。

思路分析

先把每个数 $k$ 写成标准分解形式 $\prod p_{i}^{c_{i}}$ ，注意到我们只关心 $c_{i} mod 2$ 的值，因此可以把所有的这些 $c_{i} mod 2$ 写成一个二进制数 $B_{k}$ 。

那么 $w (u, v)$ 就是 $B_{u} \oplus B_{v}$ 中所有 $1$ 对应的 $p_{i}$ 的乘积，而 $dist (u, v)$ 就是 $B_{u} \oplus B_{v}$ 中所有 $1$ 对应的 $p_{i}$ 的乘积，先逐步去掉 $B_{u} \cap (B_{u} \oplus B_{v})$ 中的 $1$ ，再逐步加入 $B_{v} \cap (B_{u} \oplus B_{v})$ 中的 $1$ ，容易证明整个过程中的数值 $\leq max (u, v)$ 。

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

容易发现 $f (p) = ⌊ \frac{n}{p} ⌋ - ⌊ \frac{n}{p^{2}} ⌋ + ⌊ \frac{n}{p^{3}} ⌋ - \dots$ ，考虑根号分治，对于 $p \leq \sqrt{n}$ 的 $p$ ，暴力计算，时间复杂度为 $O (π (\sqrt{n}) \log n) = O (\frac{\sqrt{n}}{\log n} \log n) = O (\sqrt{n})$ 。

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

因此我们只预处理在所有 $1 \sim ⌊ \frac{n}{t} ⌋$ 之间质数之和，在整除分块的过程中就可以 $O (1)$ 回答对应的区间素数数量了。

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

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

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

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

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

因此这部分的答案就是 $\sum_{p = 2}^{\sqrt{n}} φ (p) \times ⌊ \frac{n}{p^{2}} ⌋$ ，预处理出 $φ (1) \sim φ (\sqrt{n})$ 即可做到 $O (\sqrt{n})$ 。

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

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

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

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

\sum_{k = 1}^{n} μ^{2} (k) (\binom{\sqrt{⌊ n / k ⌋}}{2})

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

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

两种做法时间复杂度均为 $O (\frac{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("%lld\n",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("%lld\n",ans);
}
signed main() {
int T;
scanf("%lld",&T);
while(T--) solve();
return 0;
}