ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

[学习笔记] min_25 筛

2022-01-26 10:05:12  阅读:110  来源: 互联网

标签:25 right frac min int 质数 笔记 left mod


目录

一些约定

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

题目小结

然鹅现在只做过一道……

标签:25,right,frac,min,int,质数,笔记,left,mod
来源: https://www.cnblogs.com/AWhiteWall/p/15841383.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有