Min_25 筛笔记

给出积性函数 \(f(x)\)，要求 \(ans=\sum_{i=1}^n f(i)\)。

约定

\(lpf(i)\) 为 \(i\) 的最小质因子。

\(a/b = \lfloor \frac{a}{b} \rfloor\)

\(p_j\) 为第 \(j\) 个质数，（特殊地，\(p_0 = 1\)）

\(f'(i)\) 为质数处与 \(f(i)\) 相同的完全积性函数。

\(g(n, j) = \sum_{i=1}^n f'(i)[isprime(i) ~{\rm{or}}~ lpf(i)>p_j]\)

Part1

设 \(S(n, j) = \sum_{i=1}^n f(i)[lpf(i)>p_j]\)。

那么 \(ans = S(n, 0)+1\)

考虑拆成质数、合数部分递推得到 \(S\)：

\[\begin{aligned} S(n, j) &= g(n, |P|) - \sum_{i=1}^{j-1} f(p_i) \\ &+ \sum_{p_k^e\leq n,~k>j} f(p_k^e)(S({n}/{p_k^e}, k) + [e>1]) \end{aligned} \]

Part2

对于上面的 \(g\)，考虑递推得到：

其实 \(g(n, j)\) 就是埃筛到 \(p_j\)（第 \(j\) 轮）的时候剩下的 \([1,n]\) 没被筛掉的数 \(i\) 的 \(f'(i)\) 值贡献和。

我们有：

\[g(n, j) = g(n, j-1) ~~~~ (p_j^2>n) \\ g(n, j) = g(n, j-1) - f'(p_j)(g(n/p_j, j-1) - g(p_{j-1}, j-1)) ~~~~~ (p_j^2\leq n) \]

\(f'(p_j)(g(n/p_j, j-1) - g(p_{j-1}, j-1))\) 这部分就是 \(lpf(i)=p_j\) 的 \(i\) 的贡献。

事实上，\(g(p_{j-1}, j-1) = \sum_{i=1}^{j-1}f'(p_i)\)。

实现

• 筛出 \([1, \sqrt n ]\) 的质数以及前 \(\sqrt n\) 个 \(f\) 值
• \(f(p)\) 多项式表示的每一项筛出对应的 \(g\)，合并得到 \(O(\sqrt n)\) 个有用点值。
• 递归得到 \(S\)。

// Problem: P5325 【模板】Min_25筛
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P5325
// Memory Limit: 250 MB
// Time Limit: 2000 ms
// 
// Powered by CP Editor (https://cpeditor.org)

#include<bits/stdc++.h>
using namespace std;

#define rep(i,a,b) for(int i=(a);i<=(b);i++)
#define dwn(i,a,b) for(int i=(a);i>=(b);i--)
#define int long long

const int N=1e6+5, mod=1e9+7;

int n;
int lim;

bool vis[N];
int p[N], pc;

int s1[N], s2[N];

int add(int x, int y){
	return x+y>=mod? x+y-mod: x+y;
}

int sub(int x, int y){
	return x-y>=0? x-y: x-y+mod;
}

int mul(int x, int y){
	return x*y%mod;
}

int fpow(int x, int p){
	int res=1;
	for(; p; p>>=1, x=1LL*x*x%mod) if(p&1) res=1LL*res*x%mod;
	return res;
}

int f(int p, int c){
	// TODO: calc f(p^c)
	return mul(fpow(p, c), sub(fpow(p, c), 1));
}

const int inv2=fpow(2, mod-2), inv6=fpow(6, mod-2);

void init(){
	rep(i, 2, lim){
		if(!vis[i]){
			p[++pc]=i;
			// TODO: calc prefix sum of prime pos
			s1[pc]=add(s1[pc-1], i);
			s2[pc]=add(s2[pc-1], mul(i, i));
		}
		for(int j=1; p[j]*i<=lim; j++){
			vis[p[j]*i]=true;
			if(i%p[j]==0) break;
		}
	}
}

///////////////////////////////////////////

int val[N], id1[N], id2[N], bc;

int SUM1(int x){
	return mul(mul(x, x+1), inv2);
}

int SUM2(int x){
	return mul(mul(x, mul(x+1, x<<1|1)), inv6);
}

int g1[N], g2[N];

void get_g(){
	for(int l=1, r; l<=n; l=r+1){
		r=min(n, n/(n/l));
		val[++bc]=n/l%mod;
		
		// TODO: calc g_0
		g1[bc]=sub(SUM1(val[bc]), 1);
		g2[bc]=sub(SUM2(val[bc]), 1);
		
		val[bc]=n/l;
		if(val[bc]<=lim) id1[val[bc]]=bc;
		else id2[n/val[bc]]=bc;
	}
	
	rep(j, 1, pc){
		for(int i=1; i<=bc && p[j]*p[j]<=val[i]; i++){
			int tmp=val[i]/p[j];
			int x=(tmp<=lim? id1[tmp]: id2[n/tmp]);
			// TODO: calc g
			g1[i]=sub(g1[i], mul(p[j], sub(g1[x], s1[j-1])));
			g2[i]=sub(g2[i], mul(mul(p[j], p[j]), sub(g2[x], s2[j-1])));
		}
	}
}

int S(int i, int j){
	if(p[j]>=i) return 0;
	int x=(i<=lim? id1[i]: id2[n/i]);
	// TODO: calc Prime part
	int res=sub(sub(g2[x], g1[x]), sub(s2[j], s1[j]));
	for(int k=j+1; p[k]*p[k]<=i && k<=pc; k++){
		int pe=p[k];
		for(int e=1; pe<=i; e++, pe=pe*p[k]){
			res=add(res, mul(f(p[k], e), add(S(i/pe, k), (e>1))));
		}
	}
	return res;
}

signed main(){
	cin>>n;
	lim=sqrt(n+0.5);
	init();
	get_g();
	
	cout<<add(S(n, 0), 1);
	
	return 0;
}