Min 25筛入门

最新推荐文章于 2024-07-14 21:17:23 发布

jokerwyt

最新推荐文章于 2024-07-14 21:17:23 发布

阅读量6.2k

点赞数 7

分类专栏：新内容数论

本文链接：https://blog.csdn.net/jokerwyt/article/details/82150472

版权

新内容同时被 2 个专栏收录

86 篇文章 0 订阅

订阅专栏

数论

31 篇文章 0 订阅

订阅专栏

内容

好像是一种比较新的筛法，网上资料都是18年的
~~正式点好像就叫基于质因数分解的亚线性函数前缀和求法。~~
2020.4UPD：min25筛应用

用途

筛积性函数F的前缀和。（事实上不积性的也可以筛，只要是可以按质因子拆贡献的。）
要能构造一个完全积性函数 $G$ ，使得 $\in Prime$

复杂度

据信，空间n^0.5，时间为 $O(n^{0.75} / \ln n+n/poly(\ln n))$
（但是实际上后面那一部分还跑的挺快的）
2s大概可以做两次1e10吧。

预处理 (模拟埃氏筛法)

下文中，记 $R (a)$ 表示a的最小质因子。
一个数x要么是质数，要么 $R(x)<=\sqrt x$ ，这也是min25筛的精髓所在。他将所有数分成了两部分：质数和合数，而利用合数的上述性质，就可以快速统计答案。

首先，我们来求出所有质数的F和。
设 $g (i, j)$ 表示 $\sum F(x), 满足x\in[2,i]，且x\in Prime或R(x)>P_j$
Pj表示第j个质数。（即筛法进行到Pj时的状态）
直观地说，就是1…i做筛法筛去前j个质数后的状态。

其中j的取值只有 $\sqrt n$ 个，按j分层，我们要求的是筛法进行完毕的状态 $g (i, M A X)$
初始状态 $g [i, 0]$ 我们将所有数都视作质数来预处理，然后一个个质数去筛，将合数筛去，这样剩下的就是质数的F和。

容易推出g的递推式，用 $g (i, j - 1)$ 来推 $g (i, j)$ ，即减去所有最小质因子为 $P j$ 的数的F和。
$F(Pj)\cdot (g(i/Pj,j-1) - sumf(j-1))$
后面括号内sumf(j)是指 $\sum F(Pi),i\in[1,j]$

晒微解释一下式子的含义：
考虑任意一个需要被筛出的数 $P_j\cdot A<=i$ ，即 $A<=i / P_j$ 。
注意需要同时满足 $R(A)>=P_j$
这样的数的F和其实可以用g表示出来。但同时注意到g内有一段小于 $P j$ 的质数，需要减去。
*注意此处 $i / P j$ 可能小于Pj，此时式子不正确。但当 $Pj^2>i$ 时显然 $g (i, j) = g (i, j - 1)$ ，实现上由于省略了第二维，因此直接不更新即可。

根据常识，第一维大小只有 $2\sqrt n$ ，离散下来存好。具体可能需要看标。
下面在筛质数和。

int getid(ll x) {
	if (x <= P) return id[x]; else return id2[n/x];
}
	//p是质数数组，pre[x]表示p[2]+p[3]...+p[x]
	P = sqrt(n);
	for (ll l = 1; l <= n;) {
		ll v = n / l, r = n / v;
		if (v <= P) id[v] = ++m; else id2[n/v] = ++m; //想想为啥是单射
		w[m] = v;
		g[m] = (2+v)*(v-1)/2;
		l = r + 1;
	}
	//注意w是递减的。
	for (int j = 1; p[j] <= P; j++) {
		for (int i = 1; i <= m; i++) {
			if (p[j] * p[j] > w[i]) break; 
			g[i] -= p[j] * (g[getid(w[i] / p[j])] - pre[j-1]);
		}
	}

询问

求好了所有质数的和，就可以如何用质数和求所有和了。
其实是基于每个数的质因数分解。

有直接计算 / 递推两个方法。优缺点在最后。

直接计算版

类似g的定义，
设 $f (i, j)$ 表示 $\sum F(x), 满足x\in[2,i]，且R(x)>=P_j$
（注意是R(x)大于等于Pj而不是大于Pj）
分两个部分来计算：
1.质数
2.合数，通过枚举其最小质因子与次幂来统计。

对于f(i,j)，显然质数部分就是 $g (i, M A X) - s u m f (j - 1)$
合数部分，先枚举最小质因子与次幂 $p^k$ ，当然要满足 $p>=P_j$ 。
我们要将符合这一部分的F算进来，这样的数 $p^kA<=n且R(A)>p$ 。
由于F是积性函数，其实就是 $F(p^k)\cdot f(i/p^k,p的序号+1)$ 。递归下去求就行。出于某种神奇的原因不需要记忆化。

另外，因为f函数是不包括1的，所以质数的次幂会漏算。加上就行。

#define F(x,y) (...)
ll ask(ll x,ll i) {
	if (x <= 1 || prime[i] > x) return 0;
	ll ans = g[getid(x)] - pre[i-1];
	
	for (ll j = i; j <= prime[0] && prime[j] * prime[j] <= x; j ++) {
		ll r = prime[j];
		for (ll e = 1; r * prime[j] <= x; e++, r *= prime[j]) {
			ans = (ans + F(prime[j], e) * ask( x / r, j + 1 ) + F(prime[j], e+1)) % mo;
		}
	}
	return ans;
}

递推版

这个没有那么灵活，但是也有其优点。
使用了求g的思想。（反向模拟埃氏筛）
新设 $f (i, j)$ 表示 $\sum F(x), 满足x\in[2,i]，且x\in Prime或R(x)>=P_j$

初值 $f (i, M A X)$ 即质数和g.

$f (i, j)$ 相比于 $f (i, j + 1)$ ，需要加上
$\sum F({p_j}^kA)，满足R(A)>P_j$
同样，这样的F实际上就是 $f(i/p^k,j+1) - sumf(j)$ 再整体乘上 $F({P_j}^k)$ ，此处也漏了质数幂，需要补上。
同样在筛数的和（虽然没有什么意义…）。

	for (int j = p[0]; j; j--) {
		if (p[j] * p[j] > n) continue;
		
		for (int i = 1; i <= m; i++) {
			if (p[j] * p[j] > w[i]) break;
			ll l = p[j];
			for (int e = 1; l * p[j] <= w[i]; e++,l*=p[j]) {
				g[i] += l * (g[getid(w[i]/l)] - pre[j]) + l * p[j];
			}
		}
	}

方法比较

若只需要求一个f，那么直接计算会有比较小的常数（相比后者1/2到1/3）。
但后者可以一次处理出对于所有f(n),f(n/2),f(n/3)…f(1)的值。

步骤总结

处理出n/x的所有取值并标号，并构造一个函数G，求出对应的前缀和。
对这个G的前缀和模拟埃氏筛，求出质数的G的前缀和（事实上就是所求函数F的质数前缀和）
使用递归求法计算某个点的前缀和 / 反向模拟埃氏筛求所有n/x位置的前缀和。

例题

LOJ简单的函数
直接计算法，套式子就行。


#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
const int N = 1e5 + 1000, mo = 1e9 + 7, inv2 = 500000004;
typedef long long ll;
ll n,P;
ll id1[N],id2[N],w[2*N],m,g[2*N],h[2*N],f[2*N];
ll prime[N],pre[N];
bool is[N];

void init() {
	for (ll i = 2; i <= 1e5; i ++) {
		if (!is[i]) prime[++prime[0]] = i, pre[prime[0]] = (pre[prime[0] - 1] + i) % mo;
		for (ll j = 1; j <= prime[0] && prime[j] * i <= 100000; j ++) {
			is[prime[j] * i] = 1;
			if (i % prime[j] == 0) break;
		}
	}
}

int getid(ll x) {
	if (x <= P) return id1[x]; else return id2[n / x];
}

void calc_gh() {
	for (ll l = 1; l <= n; ) {
		ll v = n / l, r = n / v;
		if (v <= P) id1[v] = ++m; else id2[r] = ++m;
		w[m] = v;
		ll z = v % mo;
		//sum 2 ~ v
		g[m] = (2 + z) * (z - 2 + 1) % mo * inv2 % mo;
		h[m] = z - 1;
		l = r + 1;
	}
	
	for (ll i = 1; i <= prime[0]; i++) {
		ll z = prime[i];
		for (ll j = 1; j <= m && z * z <= w[j]; j++) {
			int op = getid(w[j]/z);
			g[j] = (g[j] - prime[i] * (g[op] - pre[i-1]) % mo) % mo;
			h[j] = (h[j] - (h[op] - (i - 1))) % mo;
		}
	}
	for (int i = 1; i <= m; i++) {
		g[i] = (g[i] + mo) % mo;
		f[i] = (f[i] + mo) % mo;
//		printf("%lld %lld\n",g[i],h[i]);
		f[i] = (g[i] - h[i]) % mo;
	}
}

#define F(x,y) ((x)^(y))
ll ask(ll x,ll i) {
	if (x <= 1 || prime[i] > x) return 0;
	ll ans = f[getid(x)] - (pre[i-1] - (i - 1)); if (i == 1) ans+=2;
	
	for (ll j = i; j <= prime[0] && prime[j] * prime[j] <= x; j ++) {
		ll r = prime[j];
		for (ll e = 1; r * prime[j] <= x; e++, r *= prime[j]) {
			ans = (ans + F(prime[j], e) * ask( x / r, j + 1 ) + F(prime[j], e+1)) % mo;
		}
	}
//    printf("%lld %lld %lld\n",x,i,ans);
	return ans;
}

int main() {
	cin>>n; P = sqrt(n);
	init();
	calc_gh();
	printf("%lld\n",((ask(n,1) + 1)%mo+mo)%mo);
}

jzoj5594 最大真因数
使用了递推法，方便在途中统计答案。（我是枚举最小质因子来统计的。）
（事实上只需要part1就可以完成这个统计，我是在part2统计的。）

#include <cstdio>
#include <iostream>
#include <cstring>
#include <cmath>
using namespace std;
typedef long long ll;
const int N = 2 * 750000;
ll L,R;
ll g[N+10],f[N+10],id[N+10],id2[N+10],m,w[N+10],pre[N+10];
ll p[N+10],P;
bool bz[N+10];
void init() {
	for (int i = 2; i <= N; i++) {
		if (!bz[i]) p[++p[0]] = i,pre[p[0]] = pre[p[0]-1] + i;
		for (int j = 1; j <= p[0] && i * p[j] <= N; j++) {
			bz[i * p[j]] = 1;
		}
	}
}

ll zz;
int getid(ll x) {
	if (x <= P) return id[x]; else return id2[zz/x];
}

ll calc(ll n) {
	if (n <= 1) return 1;
	zz = n;
	P = sqrt(n);
	m = 0;
	memset(id,0,sizeof id); memset(id2,0,sizeof id2);
	memset(g,0,sizeof g); memset(f,0,sizeof f);
	for (ll l = 1; l <= n;) {
		ll v = n / l, r = n / v;
		if (v <= P) id[v] = ++m; else id2[n/v] = ++m;
		w[m] = v;
		g[m] = (2+v)*(v-1)/2;
		l = r + 1;
	}
	
	for (int j = 1; p[j] <= P; j++) {
		for (int i = 1; i <= m; i++) {
			if (p[j] * p[j] > w[i]) break; //不加这句会TLE！！！ 
			g[i] -= p[j] * (g[getid(w[i] / p[j])] - pre[j-1]);
		}
	}
	
	ll ret = 0, sup = 0;
	for (int j = p[0]; j; j--) {
		if (p[j] * p[j] > n) continue;
		for (ll u = p[j]; u*p[j] <= n; u*=p[j]) {
			ret += u/p[j] * (g[getid(n/u)] - pre[j]) + u;
		}
		
		for (int i = 1; i <= m; i++) {
			if (p[j] * p[j] > w[i]) break; //不加这句会TLE！！！ 
			ll l = p[j];
			for (int e = 1; l * p[j] <= w[i]; e++,l*=p[j]) {
				g[i] += l * (g[getid(w[i]/l)] - pre[j]) + l * p[j];
			}
		}
	}
	return ret + 1;
}

int main() {
	freopen("factor.in","r",stdin);
//	freopen("factor.out","w",stdout);
	cin>>L>>R;
	init();
	cout<<calc(R) - calc(L-1)<<endl;
}

复杂度口胡

大概是 $\sum \sqrt {n/i}，并满足i<=\sqrt n$
事实上并不是每一个i都是质数，所以这是个非常宽的上界。
这个复杂度是 $\sqrt n \cdot \sum \sqrt {1/i}$ ，后面那一部分是 $n^{1/4}$ 。
证法是找一个可以求导之后等于 $i^{-0.5}$ 的函数，转换为两个端点函数值相减。
不难发现这个函数就是 $2i^{0.5}$ 。因此等于 $2(\sqrt n ^ {0.5})$ ，也就是 $O(n^{1/4})$
*为什么我觉得这个证法的界也超松的…
至于那个log怎么来的，大概是素数的缘故吧…