Min_25 筛笔记

给出积性函数 f(x),要求 ans=i=1nf(i)

约定

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

a/b=ab

pj 为第 j 个质数,(特殊地,p0=1

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

g(n,j)=i=1nf(i)[isprime(i) or lpf(i)>pj]

Part1

S(n,j)=i=1nf(i)[lpf(i)>pj]

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

考虑拆成质数、合数部分递推得到 S

S(n,j)=g(n,|P|)i=1j1f(pi)+pken, k>jf(pke)(S(n/pke,k)+[e>1])

Part2

对于上面的 g,考虑递推得到:

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

我们有:

g(n,j)=g(n,j1)    (pj2>n)g(n,j)=g(n,j1)f(pj)(g(n/pj,j1)g(pj1,j1))     (pj2n)

f(pj)(g(n/pj,j1)g(pj1,j1)) 这部分就是 lpf(i)=pji 的贡献。

事实上,g(pj1,j1)=i=1j1f(pi)

实现

  • 筛出 [1,n] 的质数以及前 nf
  • f(p) 多项式表示的每一项筛出对应的 g,合并得到 O(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;
}
posted @   HinanawiTenshi  阅读(55)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 单线程的Redis速度为什么快?
· 展开说说关于C#中ORM框架的用法!
· Pantheons:用 TypeScript 打造主流大模型对话的一站式集成库
· SQL Server 2025 AI相关能力初探
· 为什么 退出登录 或 修改密码 无法使 token 失效
点击右上角即可分享
微信分享提示