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;
}