话说我们现在要求一个函数(f)的前缀和。即求(F(n)=sum_{i=1}^nf(i))。
min25筛这个算法的主要思想是把(1...n)这些数按质数和合数分类，然后分别考虑质数和合数的贡献。
STEP1 质数贡献
我们尝试先解决一个小问题：求(G(m)=sum_{i=1}^m[iin prime]f(i))，即(m)以下的质数的(f)值之和。其中(min{frac{n}{1},frac{n}{2},...,frac{n}{n}})，共有(O(sqrt{n}))种取值。
这玩意很难直接求，于是考虑DP。设(g(j,m)=sum_{i=1}^m[iin prime ; or ; minp(i)>p[j]]f'(i))，其中(minp(i))表示(i)的最小质因子。这个式子的含义是“(i)为质数或(i)的最小质因子(>p[j])的函数值(f'(i))”之和。其中(f'(i))并不一定是我们要求的(f(i))，它可以是我们构造出来的另一个积性函数，但它必须满足以下条件：
- (f')在质数处的取值与(f)相同。即(forall pin prime)有(f'(p)=f(p))。
- (f')是完全积性函数。
- (f')可以快速求前缀和。
显然DP的边界是(g(0,m)=sum_{i=2}^mf'(i))。
考虑转移，求(g(j,m))。当(p[j]^2>m)时，(m)以内不可能有任何最小质因子是(p[j])的合数，即(p[j])无法筛掉任何数，因此(g(j,m)=g(j-1,m)quad (p[j]^2>m))。
否则，我们只需要在(g(j-1,m))里“筛掉”(leq m)的最小质因子为(p[j])的合数即可。而它们都可以表示为“(p[j] imes)一个(lfloorfrac{m}{p[j]} floor)以内的最小质因子大于等于(p[j])的数”。于是有：
这时(f')作为一个完全积性函数的好处就体现出来了，它可以直接乘以后面的一坨东西而不用担心是否互质。请读者认真理解这个式子的含义。
令(|P|)表示(leq m)的质数个数，则(g(|P|,m))就是我们前面要求的(G(m))的值了。显然，我们可以把DP的第一维滚掉。
在具体实现时，我们先预处理出(sqrt{n})以内的质数（线性筛），以及所有(m)的值（用数论分块）。暴力的做法我们可以用STL-map来记录所有(m)的值的编号，但这里可以更巧妙一些。我们开两个数组id1和id2，大小都是(sqrt{n})，巧妙地利用“对所有(>sqrt{n})的数(x)，(n/x)一定(<sqrt{n})”，就可以记下所有(m)的编号了。
int n,val[N*2],id1[N],id2[N];
//主函数（预处理）
int sqrt_n=sqrt(n),tot=0;
for(int i=1,j;i<=n;i=j+1) {
j=n/(n/i);
int w=n/i;val[++tot]=w;
if(w<=sqrt_n) id1[w]=tot;
else id2[n/w]=tot;
}
//查询某个m的编号
inline int get_id(int m) {
if(m<=sqrt_n) return id1[m];
else return id2[n/m];
}
这一部分的时间复杂度被证明是(O(frac{n^{frac{3}{4}}}{log n}))的。
STEP2 贡献加和
现在我们用刚刚搞出来的(G(m)=sum_{i=1}^m[iin prime]f(i))，求出所有(f)的和。
设(S(n,j)=sum_{i=1}^n[minp(i)geq p[j]]f(i))。则显然最终答案是(f(1)+S(n,1))。
正如一开始所说的，我们分别考虑质数和合数对(S(n,j))的贡献。显然质数的贡献是(G(n)-sum_{i=1}^{j-1}f(p[i]))。对于合数的贡献，我们可以枚举这个合数的“最小质因子及其次数”，然后进行递推。容易列出递推式：
式子里的二重循环就是在枚举合数的“最小质因子及其次数”，注意合数可以有多个不同的质因子（继续递归）或只有一个质因子的若干次方，这些都在式子里得以体现，请读者仔细理解。
在求(S)时我们直接递归，不用记忆化。这一部分的时间复杂度，被证明也是(O(frac{n^{frac{3}{4}}}{log n}))。反正你只要知道它能跑很大的数据（(1e10)的数据大约1s）就可以了。
一些栗子
预处理(f_1(i)=i^2)和(f_2(i)=i)在质数处的前缀和，在需要用到时减一下即可。
参考代码：
//P5325
#include <bits/stdc++.h>
#define ll long long
#define pb push_back
#define mk make_pair
#define pii pair<int,int>
#define fst first
#define scd second
using namespace std;
inline ll read() {
ll f=1,x=0;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}
return x*f;
}
inline void write(ll x) {
if(!x)putchar('0');if(x<0)x=-x,putchar('-');
static int sta[20];register int tot=0;
while(x)sta[tot++]=x%10,x/=10;
while(tot)putchar(sta[--tot]+48);
}
const int MAXN=1e6+5;
const int MOD=1e9+7,INV6=166666668;
ll n,sqrtn,prm[MAXN],sump[MAXN],sum2p[MAXN],cnt,g1[MAXN],g2[MAXN];
bool v[MAXN];
void sieve(int n) {
cnt=0;
for(int i=2;i<=n;++i) {
if(!v[i]) {
prm[++cnt]=i;
sump[cnt]=(sump[cnt-1]+i)%MOD;
sum2p[cnt]=(sum2p[cnt-1]+(ll)i*i%MOD)%MOD;
}
for(int j=1;j<=cnt && prm[j]*i<=n;++j) {
v[prm[j]*i]=1;
if(i%prm[j]==0) break;
}
}
}
ll val[MAXN];
int id1[MAXN],id2[MAXN];
inline int get_id(ll x){return ((x<=sqrtn)?id1[x]:id2[n/x]);}
ll getS(ll x,ll y) {
if(x<prm[y] || x<=1) return 0;
int k=get_id(x);
ll res=((g2[k]-g1[k]+MOD)%MOD-(sum2p[y-1]-sump[y-1]+MOD)%MOD+MOD)%MOD;
for(int i=y;i<=cnt && prm[i]*prm[i]<=x;++i) {
ll t1=prm[i],t2=prm[i]*prm[i];
for(int j=1;t2<=x;++j,t1=t2,t2*=prm[i]) {
ll tt1=t1%MOD,tt2=t2%MOD;
res=(res+getS(x/t1,i+1)*tt1%MOD*(tt1-1)%MOD+tt2*(tt2-1)%MOD)%MOD;
}
}
return res%MOD;
}
int main() {
ios::sync_with_stdio(0);/*syn加速*/
n=read();
sieve(sqrtn=sqrt(n));
int tot=0;/*O(sqrt(n))级别*/
for(ll i=1,j;i<=n;i=j+1) {
j=n/(n/i);
ll w=n/i;
val[++tot]=w;
if(w<=sqrtn) id1[w]=tot;
else id2[n/w]=tot;
w%=MOD;
g1[tot]=w*(w+1)/2%MOD;
g2[tot]=w*(w+1)%MOD*(2LL*w+1)%MOD*INV6%MOD;
g1[tot]=(g1[tot]-1+MOD)%MOD;
g2[tot]=(g2[tot]-1+MOD)%MOD;
}
for(int j=1;j<=cnt;++j) {
for(int i=1;i<=tot && prm[j]*prm[j]<=val[i];++i) {
int k=get_id(val[i]/prm[j]);
g1[i]=(g1[i]-prm[j]*(g1[k]-sump[j-1]+MOD)%MOD+MOD)%MOD;
g2[i]=(g2[i]-prm[j]*prm[j]%MOD*(g2[k]-sum2p[j-1]+MOD)%MOD+MOD)%MOD;
}
}
//for(int i=1;i<=tot;++i)cout<<val[i]<<" "<<g1[i]<<" "<<g2[i]<<endl;
write((getS(n,1)+1)%MOD);puts("");
return 0;
}
先看(sum_{i=1}^nvarphi(i))。首先(varphi(p)=p-1quad (pin prime))。于是我们预处理(f_1(i)=i)和(f_2(i)=1)的前缀和，按与上一题相同的方法相减即可求出(G(m))。
另外，(varphi(p^c)=p^{c-1}(p-1))（自己yy可得），这样在求(S)时，质数的整次幂的函数值也很好求。
然后是(mu(i))，显然(mu(p)=-1quad (pin prime))。于是我们只要对前面已经求好的(f_2(i)=1)求个相反数即可。
另外，(mu(p^k)=0quad (k>1))，因此求(S)时不需要枚举最小质因子的次数。
参考代码：
//P4213 min_25
#include <bits/stdc++.h>
#define ll long long
#define pb push_back
#define mk make_pair
#define pii pair<int,int>
#define fst first
#define scd second
using namespace std;
/* ----- by:duyi ----- */
const int N=5e4;
const int MAXN=N*2+5;
int n,sn,cnt,tot,p[MAXN],sum[MAXN],id1[MAXN],id2[MAXN],val[MAXN];
ll g1[MAXN],g2[MAXN];
bool v[MAXN];
void sieve() {
v[1]=1;
for(int i=2;i<=N;++i) {
if(!v[i]) p[++cnt]=i;
for(int j=1;j<=cnt && (ll)i*p[j]<=N;++j) {
v[i*p[j]]=1;
if(i%p[j]==0) break;
}
}
for(int i=1;i<=cnt;++i) sum[i]=sum[i-1]+p[i];
}
inline int get_id(int x) {
if(x<=sn) return id1[x];
else return id2[n/x];
}
ll S_phi(int x,int y) {
if(x<=1 || p[y]>x) return 0;
ll res=g1[get_id(x)]-g2[get_id(x)]-(sum[y-1]-(y-1));
for(int i=y;i<=cnt && (ll)p[i]*p[i]<=x;++i) {
ll pre=1,cur=p[i];
for(int j=1;cur*p[i]<=x;++j) {
res+=pre*(p[i]-1LL)*S_phi(x/cur,i+1)+cur*(p[i]-1LL);
pre=cur;cur=cur*p[i];
}
}
return res;
}
ll S_mu(int x,int y) {
if(x<=1 || p[y]>x) return 0;
ll res=-g2[get_id(x)]+y-1;
for(int i=y;i<=cnt && (ll)p[i]*p[i]<=x;++i) {
res+=(-S_mu(x/p[i],i+1));
}
return res;
}
void solve() {
cin>>n;
sn=sqrt(n);tot=0;
for(int i=1,j;i<=n;i=j+1) {
j=n/(n/i);int w=n/i;
val[++tot]=w;
if(w<=sn) id1[w]=tot;
else id2[n/w]=tot;
g1[tot]=(ll)w*(w+1LL)/2LL-1LL;
g2[tot]=w-1;
}
for(int i=1;i<=cnt;++i) {
for(int j=1;j<=tot && (ll)p[i]*p[i]<=val[j];++j) {
int t=get_id(val[j]/p[i]);
g1[j]-=(ll)p[i]*(g1[t]-sum[i-1]);
g2[j]-=(g2[t]-(i-1));
}
}
cout<<S_phi(n,1)+1LL<<" "<<S_mu(n,1)+1LL<<endl;
}
int main() {
ios::sync_with_stdio(0);/*syn加速*/
sieve();
int T;cin>>T;while(T--)solve();
return 0;
}
总结一下
用min_25做题主要是要想清楚两个问题：
-
该函数在质数处的取值怎么求和？（找到合适的(f')）
-
质数的次幂处的取值(f(p^c))能否快速求？
本文写得比较粗浅，欢迎大家补充。如有纰漏欢迎指正。