设 $f$ 为积性函数,且其在质数处的取值 $f(p)$ 是关于 $p$ 的多项式,且其在质数的幂的取值 $f(p^c)$ 能快速计算,则 Min_25 筛可以在 $\mathcal{O}(\frac{n^{\frac{3}{4}}}{\log n})$ 的时间复杂度内算出 $\sum_{i=1}^nf(i)$。
思路
为了方便,设 $\mathrm{P}$ 为质数集,$p_i$ 为第 $i$ 小的质数。
第一步
先考虑如何求 $\sum_{i=1}^n[i\in\mathrm{P}]f(i)$。
因为 $f(p)$ 是关于 $p$ 的多项式,因此我们考虑对于若干个函数 $f_k(i)=i^k$,计算 $\sum_{i=1}^x[i\in\mathrm{P}]f_k(i)$,然后再加减一下就可以得到 $\sum_{i=1}^x[i\in\mathrm{P}]f(i)$ 了。
考虑怎么求 $\sum_{i=1}^x[i\in\mathrm{P}]f_k(i)$。为了方便,本节剩余部分中直接将 $f_k(i)$ 写作 $f(i)$。
考虑一个 DP。设
$$
g(n,j)=\sum_{i=1}^n[i\in\mathrm{P}\lor\min_{p|i,p\in\mathrm{P}}\left\{p\right\}>p_j]f(i)
$$
可以发现,$\sum_{i=1}^n[i\in\mathrm{P}]f(i)$ 就等于 $g(n,|\mathrm{P}|)$。
考虑 $g(n,j)$ 的实际意义。可以发现 $g(n,j)$ 就等于埃筛时用 $j$ 个质数筛后未被筛去的所有数的 $f$ 值之和加上所有质数的 $f$ 值之和。
根据这个实际意义考虑转移。
- 当 $p_j^2>n$ 时,显然不会筛去任何数,所以此时 $g(n,j)=g(n,j-1)$。
- 当 $p_j^2\leq n$ 时,筛掉了最小质因子为 $p_j$ 的那些数。考虑将这些数除掉一个 $p_j$,那么除完之后的数最小质因子要 $\geq p_j$,则这些数的 $f$ 值之和为 $g(\frac{n}{p_j},j-1)-\sum_{i=1}^{j-1}f(p_i)$。又因为 $f$ 是完全积性函数,所以此时 $g(n,j)=g(n,j-1)-f(p_j)\times\left[g(\frac{n}{p_j},j-1)-\sum_{i=1}^{j-1}f(p_i)\right]$。
综上可以得到总转移
$$
g(n,j)=\begin{cases}g(n,j-1),&p_j^2>n\\g(n,j-1)-f(p_j)\times\left[g(\frac{n}{p_j},j-1)-\sum_{i=1}^{j-1}f(p_i)\right],&p_j^2\leq n\end{cases}
$$
边界为 $g(n,0)=\sum_{i=1}^nf(i)$。
第二步
设
$$
S(n,j)=\sum_{i=1}^n\left[\min_{p|i,p\in\mathrm{P}}\left\{p\right\}\geq p_j\right]f(i)
$$
考虑转移,分质数和合数两部分计算,可以得到
$$
S(n,j)=g(n,|\mathrm{P}|)-\sum_{i=1}^{j-1}f(p_i)+\sum_{k\geq j,p_k^2\leq n}\sum_{e\geq 1,p_k^{e+1}\leq n}\left[S\left(\frac{n}{p_k^e},k+1\right)\times f(p_k^e)+f(p_k^{e+1})\right]
$$
合数部分的意思是,枚举最小质因数以及次数,由于是积性函数所以可以直接乘,然后少算了 $f(p_k^2),f(p_k^3),\cdots$ 所以需要补进来($f(p_k)$ 在质数部分计算过了)。
那么最后答案即为 $S(n,1)+f(1)$。
细节
递归求解 $S(n,j)$。则递归到的那些 $n$ 为 $\left\lfloor\frac{n}{a}\right\rfloor,\left\lfloor\frac{\left\lfloor\frac{n}{a}\right\rfloor}{b}\right\rfloor,\cdots$。然后有一个这样的定理
$$
\left\lfloor\frac{\left\lfloor\frac{n}{a}\right\rfloor}{b}\right\rfloor=\left\lfloor\frac{n}{ab}\right\rfloor
$$
因此递归到的 $n$ 为 $\left\lfloor\frac{n}{i}\right\rfloor$,而这样的值的数量是 $\mathcal{O}(\sqrt{n})$ 级别的。
于是我们仅需要计算出 $\mathcal{O}(\sqrt{n})$ 个 $g(n,|\mathrm{P}|)$ 即可。
计算 $g(n,|\mathrm{P}|)$ 时,当 $p_j>n$ 时 $g$ 值不会再改变,因此只需要预处理出 $\sqrt{n}$ 内的素数即可。
实现
首先筛出 $\sqrt{n}$ 内的质数。
然后求出所有 $\left\lfloor\frac{n}{i}\right\rfloor$ 的值并存下来。存的时候开两个数组,如果 $\left\lfloor\frac{n}{i}\right\rfloor\leq\sqrt{n}$ 则以 $\left\lfloor\frac{n}{i}\right\rfloor$ 为下标存到第一个数组中,否则以 $\left\lfloor\frac{n}{\left\lfloor\frac{n}{i}\right\rfloor}\right\rfloor$ 为下标存到第二个数组中。
接着对求出的这些值,求出它对应的 $g(n,|\mathrm{P}|)$。
最后直接递归求解 $S(n,1)$ 即可,记得加上 $f(1)$。
理解
本部分不一定准确,如有错误烦请指正。
由于最后算出来的 $g(n,|\mathrm{P}|)$ 仅包含质数处的取值,所以我们假装 $f$ 在非质数处的取值和质数处相同是没有影响的。
在计算 $g(n,|\mathrm{P}|)$ 时的转移要求 $f$ 是完全积性函数,所以需要拆成若干个函数分开求。这就是要求 $f(p)$ 为关于 $p$ 的多项式的原因。
计算 $S(n,j)$ 时的转移要求 $f$ 是积性函数且 $f(p^c)$ 能够快速计算。
例子
$f(p)=p^2-p$,因此分别对 $f_1(i)=i$ 和 $f_2(i)=i^2$ 求 $g(n,|\mathrm{P}|)$,答案为 $S(n,1)+1$。
// ===================================
// author: M_sea
// website: http://m-sea-blog.com/
// ===================================
#include <bits/stdc++.h>
#define re register
using namespace std;
typedef long long ll;
inline ll read() {
ll X=0,w=1; char c=getchar();
while (c<'0'||c>'9') { if (c=='-') w=-1; c=getchar(); }
while (c>='0'&&c<='9') X=X*10+c-'0',c=getchar();
return X*w;
}
const int N=200000+10;
const int mod=1e9+7,inv6=166666668,inv2=500000004;
ll n; int s;
int p[N],v[N],cnt=0;
inline void sieve(int n) {
for (re int i=2;i<=n;++i) {
if (!v[i]) p[++cnt]=i;
for (re int j=1;j<=cnt&&i*p[j]<=n;++j) {
v[i*p[j]]=1;
if (i%p[j]==0) break;
}
}
}
ll w[N]; int id1[N],id2[N]; int top=0;
int gs1[N],gs2[N],g1[N],g2[N];
inline int S(ll x,int j) {
if (x<=1||p[j]>x) return 0;
int k=x<=s?id1[x]:id2[n/x];
int res=((1ll*g2[k]-g1[k]-gs2[j-1]+gs1[j-1])%mod+mod)%mod;
for (re int i=j;i<=cnt&&1ll*p[i]*p[i]<=x;++i) {
ll p1=p[i],p2=1ll*p[i]*p[i];
for (re int e=1;p2<=x;++e,p1=p2,p2*=p[i]) {
int q1=p1%mod,q2=p2%mod;
res=(res+1ll*S(x/p1,i+1)*q1%mod*(q1-1)+1ll*q2*(q2-1))%mod;
}
}
return res;
}
int main() {
n=read(),s=sqrt(n); sieve(s);
for (re int i=1;i<=cnt;++i) {
gs1[i]=(gs1[i-1]+p[i])%mod;
gs2[i]=(gs2[i-1]+1ll*p[i]*p[i])%mod;
}
for (re ll l=1,r;l<=n;l=r+1) {
r=n/(n/l); w[++top]=n/l;
n/l<=s?id1[n/l]=top:id2[n/(n/l)]=top;
int x=(n/l)%mod;
g1[top]=(1ll*x*(x+1)%mod*inv2%mod-1+mod)%mod;
g2[top]=(1ll*x*(x+1)%mod*(x+x+1)%mod*inv6%mod-1+mod)%mod;
}
for (re int i=1;i<=cnt;++i)
for (re int j=1;j<=top&&1ll*p[i]*p[i]<=w[j];++j) {
ll x=w[j]/p[i]; int k=x<=s?id1[x]:id2[n/x];
g1[j]=(g1[j]-1ll*p[i]*(g1[k]-gs1[i-1]+mod)%mod+mod)%mod;
g2[j]=(g2[j]-1ll*p[i]*p[i]%mod*(g2[k]-gs2[i-1]+mod)%mod+mod)%mod;
}
printf("%d\n",(S(n,1)+1)%mod);
return 0;
}