\[G(n,k)=\begin{cases}G(n,k-1)&p_k^2>n\\G(n,k-1)-f"(p_k)\times\left(G(\lfloor n/p_k\rfloor,k-1)-\sum_{i=1}^{k-1}f"(p_i) \right) &p_k^2\le n\end{cases}\]同样在递归的过程中也只会访问到所有的 \(\lfloor n/t\rfloor\) 下标处，且最后 \(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/t\rfloor )\) 处的值，应该怎么处理呢？
最简单的想法是在递归求 \(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)\) 表示 \(1\sim n\) 中，满足 \(i\) 的最小质因子 \(>p_k\) 或 \(i\) 是质数的 \(f(i)\) 的和。
依然考虑 \(S(n,k-1)-S(n,k)\)，在 \(p_k^2> n\) 时依然为 \(0\)，在 \(p_k^2\le n\) 时枚举 \(p_k\) 的次数得到：
\[S(n,k-1)-S(n,k)=\sum_{c=1}^{p_k^c\le n} f(p_k^c)\times\left(S(\lfloor n/p_k^c\rfloor,k)+[c>1]-\sum_{i=1}^{p_i\le \min(p_k,\lfloor n/p_k^c\rfloor)}f(p_i)\right)\]但是注意到最后一项求和式并不是很好看，考虑化简，注意到当 \(p_k> \lfloor n/p_k^c\rfloor\) 时，\(p_k^{c+1}> n\)，那么 \(S(\lfloor n/p_k^c\rfloor,k)\) 只会统计 \(\le \lfloor n/p_k^c\rfloor\) 的质数，与求和式正好抵消，因此可以把上式写成：
\[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)\times\left(S(\lfloor n/p_k^c\rfloor,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)+\sum\limits_{c=1}^{p_k^{c+1}\le n}f(p_k^{c+1})+ f(p_k^c)\times\left(S(\lfloor n/p_k^c\rfloor,k)-\sum_{i=1}^{k}f(p_i)\right)&p_k^2\le n\end{cases}\]从大到小枚举 \(p_k\) 再从大到小枚举 \(n\) 转移即可。
其中边界条件为 \(S(n,\infty)=Q(n)\)，最终的答案是 \(S(\lfloor n/t\rfloor,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/i\rfloor)\) 的和式计算问题，其中 \(f(i)\) 是积性函数，\(g(i)\) 是任意一个关于 \(i\) 的函数。
那么我们根据整除分块枚举 $\lfloor n/i\rfloor $，则 \(g\) 函数的值固定，要求的就是某一段特定区间内 \(f(i)\) 的和，容易发现这些区间的端点都是某些 \(\lfloor n/t\rfloor\)，用上面所介绍的非递归型的 Min-25 筛处理即可。
三、代码实现以线性筛模板题为例（Link）下面的代码是递归的实现：
#include#define int long longusing 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#define int long longusing 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;}
实际运行上，非递归版本的实现会慢于递归版的实现。
四、复杂度分析空间复杂度 \(\mathcal O(\sqrt n)\)，可以参照参考代码的实现。
时间复杂度可以看做合法的 \(S(n,k) / G(n,k)\) 状态数，对于每个 \(n=\lfloor n/t\rfloor\)，合法的 \(k\) 共有 \(\pi(\sqrt{\lfloor n/t\rfloor})\) 个。
因此复杂度可以估计为：
\[\begin{aligned}\mathcal T(n)&=\sum_{i^2\le n} \mathcal O(\pi(\sqrt i))+\mathcal O(\pi(\sqrt{n/i}))\\[2.5ex]&=\sum_{i^2\le n} \mathcal O\left(\dfrac{\sqrt i}{\ln \sqrt i}\right)+\mathcal O\left(\dfrac{\sqrt{n/i}}{\ln \sqrt{n/i}}\right)\\[2.5ex]&=\mathcal O\left(\dfrac{1}{\log n}\int_1^{\sqrt n} \sqrt{\dfrac nx} \mathrm dx\right)\\[2.5ex]&=\mathcal O\left(\dfrac{\sqrt n}{\log n}\times\left. 2\sqrt x \right|^{\sqrt n}_1 \right)\\[2.5ex]&=\mathcal O(\dfrac{n^{3/4}}{\log n})\end{aligned}\]五、例题I. [SPOJ-34096] DIVCNTKProblem Link
题目大意
给定 \(n\)，求 \(\sum_{i=1}^n \sigma_0(i^k)\)（即 \(i^k\) 的因子个数）。
数据范围：\(n\le 10^{10}\)。
思路分析
容易证明 \(\sigma_0(i)\) 是积性函数，同理 \(\sigma_0(i^k)\) 同样也是积性函数，考虑求 \(\sigma_0(p^k)\)，显然 \(\sigma_0(p^k)=k+1\)，我们只需要求所有 \(1\sim \lfloor n/t\rfloor\) 之间的质数个数，构造 \(f"(p)=1\)，直接用 Min-25 筛处理即可。
时间复杂度 \(\mathcal O(\dfrac{n^{3/4}}{\log n})\)
代码呈现
#include#define int unsigned long longusing 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
II. [洛谷-4213] SUMProblem Link
题目大意
求 \(\mu(i),\varphi(i)\) 的前缀和。
数据范围：\(n\le 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#define int long longusing 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 APSPProblem Link
题目大意
给定一张 \(n\) 个点的无向完全图，\(i,j\) 之间的边权定义为最小的正整数 \(x\) 使得 \(ijx\) 是完全平方数。
求 \(\sum_{1\le i\le j\le n}\mathrm{dist}(i,j)\) 的值，其中 \(\mathrm{dist}(i,j)\) 表示 \((i,j)\) 之间的最短路长度。
数据范围：\(n\le 10^{10}\)。
思路分析
先把每个数 \(k\) 写成标准分解形式 \(\prod p_i^{c_i}\)，注意到我们只关心 \(c_i\bmod 2\) 的值，因此可以把所有的这些 \(c_i\bmod 2\) 写成一个二进制数 \(B_k\)。
那么 \(w(u,v)\) 就是 \(B_u\oplus B_v\) 中所有 \(1\) 对应的 \(p_i\) 的乘积，而 \(\mathrm{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\)，容易证明整个过程中的数值 \(\le \max(u,v)\)。
因此考虑拆贡献，对于某个质数 \(p\)，记 \(f(p)\) 为 \([1,n]\) 中含有奇数个 \(p\) 的数的数量，那么 \(p\) 对答案的贡献就是 \(p\times f(p)\times(n-f(p))\)。
容易发现 \(f(p)=\left\lfloor\dfrac{n}{p}\right\rfloor-\left\lfloor\dfrac{n}{p^2}\right\rfloor+\left\lfloor\dfrac{n}{p^3}\right\rfloor-\dots\)，考虑根号分治，对于 \(p\le \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)=\left\lfloor\dfrac np\right\rfloor\)，用整除分块的技巧枚举 \(\left\lfloor\dfrac np\right\rfloor=k\)，那么 \(f(p)\times (n-f(p))\) 容易计算，那么我们只要求某个特定区间 \([l,r]\) 之中的质数之和，考虑整除分块的过程，注意到 \(l-1,r\) 都能写成某个 \(\left\lfloor\dfrac nt\right\rfloor\) 的形式。
因此我们只预处理在所有 \(1\sim \left\lfloor\dfrac nt\right\rfloor\) 之间质数之和，在整除分块的过程中就可以 \(\mathcal O(1)\) 回答对应的区间素数数量了。
容易发现这是一个标准的 Min-25 筛的第一部分，因此直接用 Min-25 筛处理即可。
这一部分的时间复杂度为 \(\mathcal O(\dfrac{n^{3/4}}{\log n})\)。
但是我们还漏掉了一些贡献，对于 \(u\times v\) 是完全平方数但 \(u\ne v\) 的数对 \((u,v)\) 对答案是有 \(1\) 的贡献的，因此我们要求，有多少对 \((u,v)\) 满足 \(1\le u设公比最简分数为 \(q/p\)，枚举分母 \(p\)，此时分子 \(q然后考虑 \(v\)，注意到 \(u=\dfrac {vq^2}{p^2}\)，因此 \(p^2\mid v\)，这样的 \(v\) 共有 \(\left\lfloor\dfrac n{p^2}\right\rfloor\) 种，容易证明这样的每一对 \((p,q,v)\) 都和我们要计数的 \((u,v)\) 构成双射。
因此这部分的答案就是 \(\sum_{p=2}^{\sqrt n} \varphi(p)\times \left\lfloor\dfrac n{p^2}\right\rfloor\)，预处理出 \(\varphi(1)\sim \varphi(\sqrt n)\) 即可做到 \(\mathcal O(\sqrt n)\)。
将答案分成三部分分别处理即可。
第三部分的另一种处理方法
对于这样乘积是完全平方数的数对计数，最朴素的思路就是枚举 \(u,v\) 在各自除掉最大的完全平方因子后得到的数 \(k\)，枚举这样的 \(k\)，\(u,v\) 的选择方案数就是 \(g(\lfloor n/k\rfloor)=\binom{\sqrt{\lfloor n/k\rfloor}}2\)。
考虑刻画 \(k\) 以写出求和式，显然 \(k\) 中无平方因子，则 \(\mu(k)\ne 0\)，因此满足条件的 \(k\) 可以用 \(\mu^2(k)=1\) 表示，因此这部分的贡献就是：
\[\sum_{k=1}^n \mu^2(k)\binom{\sqrt{\lfloor n/k\rfloor}}2\]考虑整除分块枚举 \(\lfloor n/k\rfloor\)，原问题变成统计 \(\mu^2(k)\) 在所有 \(\lfloor n/t\rfloor\) 处的前缀和，注意到 \(\mu^2(k)\) 是积性函数，因此可以用非递归版本的 Min-25 筛求出。
由于 \(\mu^2(k)\) 函数的特殊性，在递推 \(S(n,k)\) 的时候只需要处理 \(c=1\) 的情况，因此该算法的运行效率和前一个算法的运行效率差别并不大。
两种做法时间复杂度均为 \(\mathcal O(\dfrac{n^{3/4}}{\log n})\)。
代码呈现
实现方式一：
#include#define int long long#define LL __int128using 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
实现方式二：
#include#define int long long#define LL __int128using 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;}
关键词：
(相关资料图)