[学习笔记] min_25 筛

目录

一些约定

  1. \(p_k\) 表示第 \(k\) 个质数,特别地,\(p_0=1\);
  2. 令 \(n/k\) 为 \(\left \lfloor \frac{n}{k}\right \rfloor\);
  3. 令 \(\text{lpf}(x)\) 为数 \(x\) 的最小质因数。

要干啥?

积性函数 \(F(i)\) 的前缀和。

推柿子

令 \(h(i,j)=\sum_{x=2}^i[\text{lpf}(x)> p_j]F(x)\),那么答案就是 \(h(n,0)+F(1)=h(n,0)+1\)(积性函数的第一项一定为 \(1\),\(1\) 和任何数互质嘛)。求解的普遍思路似乎是递推,不妨先按照定义展开(注意,本次展开只包含可被分解成多个质数的合数):

\[h(i,j)=\sum_{k=j+1}\sum_{e=1} F(p_k^e)\cdot h(i/p_{k}^e,k) \]

其实就是枚举最小质因子及其幂次,然后进行递归(易得它们是互质的)。

另外就还剩下质数、单个质数的幂的情况,后者也是比较容易写的:

\[h(i,j)=\sum_{k=j+1}\sum_{e=1} F(p_k^e)\cdot h(i/p_{k}^e,k)+F(p_k^{e+1})\ \ \ \ (p_{k}^{e+1}\le i) \]

需要注意后面的限制,可以发现,当 \(p_k^{e+1}>i\) 时,\(F(p_k^e)\cdot h(i/p_{k}^e,k)\) 一定也为零,所以这个限制是紧凑且合理的。

令 \(P(n)\) 为小于等于 \(n\) 的质数的 \(F\) 值的前缀和,那么最终就有(注意到质数只需要枚举到 \(\sqrt n\) 级别):

\[h(i,j)=P(i)-P(p_j)+\sum_{k=j+1}\sum_{e=1} F(p_k^e)\cdot h(i/p_{k}^e,k)+F(p_k^{e+1})\ \ \ \ (p_{k}^{e+1}\le i) \]

于是我们只用计算 \(P(n)\).

首先,我们需要 \(F(n)\) 是一个 次数较小 的多项式,这样拆开的每一项就可以视作一个 完全积性函数 —— 更一般地,不妨设 \(f_k(n)=n^k\),只需要计算这个函数的 \(P_k(n)\),乘上系数再相加即可得到 \(P(n)\).

令 \(g(i,j)\) 为在 \([2,i]\) 之间且 \(\text{lpf}>p_j\) 的数字的 \(f\) 值之和,特别地,这一坨恒包含小于等于 \(i\) 的质数的 \(f\) 值。

考虑递推地求取 \(g(i,j)\) —— 我们可以将 \(g(i,j)\) 理解为在 \(g(i,j-1)\) 的基础上,去掉 最小质因子 为 \(p_j\) 的 合数 的 \(f\) 值:

\[g(i,j)=g(i,j-1)-f_k(p_j)\cdot\left (g(i/p_j,j-1)-\sum_{q=1}^{j-1}f_k(p_q)\right)\ \ \ \ \ (i\ge p_j^2) \]

由于这是完全积性函数,那么有 \(f(r\cdot p_j)=f(r)\cdot f(p_j)\),后面的一大坨相当于计算合法的 \(r\). 可以发现完全积性函数的性质省去了枚举质数幂次的操作。

注意到后面 \(i\ge p_j^2\) 的条件。实际上,若 \(i<p_j^2\),那么就不可能找到合法的 \(r\),所以直接 \(g(i,j)\leftarrow g(i,j-1)\) 即可。

容易发现,\(P_k(n)=g(n,cnt)\)(\(cnt\) 是筛出质数的个数)。


关于复杂度:\(\text{It is said that}\) 计算 \(h(i,j)\) 的复杂度为 \(\mathcal O(n^{1-\epsilon})\),计算 \(g(i,j)\) 的复杂度为 \(\mathcal O\left(\frac{n^{\frac{3}{4}}}{\log n}\right)\). 证明不会,可能会一直咕下去。

代码实现

参照某谷 模板题

在实现过程中,函数的 \(i,j\) 参数的 个数 都是小于 \(\sqrt n\) 的。\(j\) 的原因就不阐述,对于 \(i\),可以观察到它是一个类似 \(\left \lfloor \frac{\left \lfloor \frac{n}{a}\right \rfloor}{b}\right \rfloor\) 的向下嵌套的结构,可以证明它的取值属于 \(\left \lfloor \frac{n}{a}\right \rfloor\) 的取值集合,这样不仅可以证明它的个数在 \(\sqrt n\) 之内,而且还可以进行整除分块预处理。以下给出感性证明:不妨将 \(\left \lfloor \frac{n}{a}\right \rfloor\) 中的分母替换成 \(a'\) 使得 \(a'\mid n\) 且 \(\frac{n}{a'}=\left \lfloor \frac{n}{a}\right \rfloor\),这样嵌套就可以换成 \(\left \lfloor \frac{n}{a'b}\right \rfloor\),结论得证。

于是我们将合理的 \(i\) 存起来,需要注意的是,\(i\) 的取值可能很大,所以对于 \(i\le \sqrt n\) 的部分用 \(i\) 来查找,另一部分用 \(n/i\) 来查找(这个玩意实际上就是整除分块时有 \(i=n/l\),对应的 \(r\) 就是 \(n/i\))。

回到此题,可知 \(F(n)=n^2-n\),可以分成 \(f_1,f_2\) 来计算,然后合并即可。

#include <cstdio>
#define print(x,y) write(x),putchar(y)

template <class T>
inline T read(const T sample) {
	T x=0; char s; bool f=0;
	while((s=getchar())>'9' || s<'0')
		f |= (s=='-');
	while(s>='0' && s<='9')
		x = (x<<1)+(x<<3)+(s^48),
		s = getchar();
	return f?-x:x;
}

template <class T>
inline void write(T x) {
	static int writ[50],tp=0;
	if(x<0) putchar('-'),x=-x;
	do writ[++tp] = x-x/10*10, x/=10; while(x);
	while(tp) putchar(writ[tp--]^48);
}

#include <cmath>
#include <iostream>
using namespace std;
typedef long long ll;

const int maxn = 1e6+5;
const int mod = 1e9+7;
const int inv = 166666668;

bool is[maxn];
ll n,id[maxn];
int S,f[2][maxn],pc,p[maxn];
int g[2][maxn],tot,ind[2][maxn];

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

inline int dec(int x,int y) {
	return x-y<0?x-y+mod:x-y;
}

void sieve(int n) {
	for(int i=2;i<=n;++i) {
		if(!is[i]) {
			p[++pc] = i;
			f[0][pc] = inc(f[0][pc-1],i);
			f[1][pc] = inc(f[1][pc-1],1ll*i*i%mod);
		}
		for(int j=1; j<=pc && i*p[j]<=n; ++j) {
			is[i*p[j]] = 1;
			if(i%p[j]==0) break;
		}
	}
}

inline int getSum(ll n,bool d) {
	int ret; n %= mod; // qwq
	if(!d) ret = (n*(n+1)>>1)%mod;
	else ret = n*(n+1)%mod*(n*2+1)%mod*inv%mod;
	return (ret-1<0?ret-1+mod:ret-1);
}

void getf() {
	ll l,r,x;
	for(l=1;l<=n;l=r+1) {
		r = min(n,n/(n/l)); 
		id[++tot] = (x=n/l);
		if(x<=S) ind[0][x]=tot;
		else ind[1][n/x]=tot; // for better restoration
		g[0][tot] = getSum(x,0);
		g[1][tot] = getSum(x,1);
	}
    // 用滚动数组递推 g(i,j)
	for(int i=1;i<=pc;++i) {
		for(int j=1; j<=tot && 1ll*p[i]*p[i]<=id[j]; ++j) {
			int pre = (id[j]/p[i]<=S)?
					ind[0][id[j]/p[i]]:ind[1][n/(id[j]/p[i])];
			g[0][j] = dec(g[0][j],1ll*p[i]*dec(g[0][pre],f[0][i-1])%mod);
			g[1][j] = dec(g[1][j],1ll*p[i]*p[i]%mod*dec(g[1][pre],f[1][i-1])%mod);
		}
	}
}

int h(ll x,int y) {
	if(x<=p[y]) return 0;
	int ID = (x<=S)?ind[0][x]:ind[1][n/x];
	int ans = dec(dec(g[1][ID],g[0][ID]),dec(f[1][y],f[0][y]));
	for(int i=y+1; i<=pc && 1ll*p[i]*p[i]<=x; ++i) {
		for(ll tmp=p[i];;tmp=tmp*p[i]) {
			if(tmp*p[i]>x) break;
			ans = inc(ans,inc(
				tmp*(tmp-1)%mod*h(x/tmp,i)%mod,
				tmp*p[i]%mod*((tmp*p[i]-1)%mod)%mod
			));
		}
	}
	return ans;
}

int main() {
	n=read(9ll);
	sieve(S=sqrt(n)); getf(); 
	print(inc(h(n,0),1),'\n');
	return 0;
}

题目小结

然鹅现在只做过一道……

上一篇:python中上下文管理,with的用法


下一篇:MaskEdit 控件IP地址的Mask设置