[学习笔记] Min_25筛

这个东西 ( t jzm) 又会了!但是我自己看了一会还是看懂了。

什么是Min25?

据说是一个叫 ( t Min25) 的人发明的,跟杜教筛的玄学程度有得一拼,不过并不需要很多的前置知识。

它是用来解决这样一类问题:已知 (f(p^k)) 是关于质数 (p) 的多项式,(f(x)) 是积性函数,求 (sum_{i=1}^n f(i)) ,做法是把积性函数拆成若干个完全积性函数,算出对于 (n) 以内质数的函数值,再用递归的形式算答案,总的来说,体现了用 函数结构,找递归子问题 来加速计算的思想(再杜教筛的题里面也很常见)。

如果你以前没有学过 ( t Min25) ,你肯定听不懂我在讲什么,没关系,我们慢慢来解析这个算法。

由于 (1) 很特殊,所以我们下面的算法是不考虑 (1) 的,最后把 (1) 的答案加上即可。

分类

首先对 (i) 按照质数和合数分类:

[sum_{i=1}^nf(i)=sum_{1leq pleq n}f(p)+sum_{ t otherwise}^n f(i) ]

对于后面的数枚举 最小质因数 ,因为合数的最小质因数 (leqsqrt n) 所以会快(这也印证了为什么单独讨论质数):

[sum_{1leq pleq n}f(p)+sum_{p,e}f(p^e)Big(sum_{i=1,minp>p}^{n/p^e}f(i)Big) ]

其中 (minp) 表示 (i) 的最小质因数,发现了吗?我们把质数提出来后,利用积性函数的性质,把最小质因子拆出来,就得到了一个子问题 ,这时候我们的思路就渐渐明晰了,算法分为两部分:处理前面质数的答案,处理后面的递归子问题。

第一部分

这是 ( t Min25) 的神奇之处,我们可以在不知道所有质数的情况下求出所有质数的答案之和。

(g(n,j)) 表示 (n) 以内的数,满足是个质数或者最小质因数 (>p_j) 的合数的函数值,可能我们需要的函数值是一个多项式,但是这里我们把他改成一个形如 (x^k) 的东西,这样他就是一个 完全积性函数 ,因为 (g(..,0)) 是好算的(直接求和),然后我们用 (sqrt n) 以内的质数慢慢 把合数筛掉 ,最后就只剩质数了,我们来看一看怎么筛的:

[g(n,j)=g(n,j-1)-p_j^k(g(frac{n}{p_j},j-1)-sp(j-1)) ]

减掉的就是包含这个质因子,但是不是最小质因子的数,因为是完全积性函数所以不用把所有的质因子提出来,而只用提一个。后面减掉的是 (p_{j-1}) 以内的质数,这些是不应该算进去的,就是 (sp) 数组,它可以线性筛出来。

第一维很大怎么办,但是你发现我们用到的只有 (frac{n}{x}) 这些取值,一共只有 (sqrt n) 个,所以会很快。

第二部分

要解决递归的问题了,设 (S(n,j)) 表示 (n) 以内的数最小质因子大于 (p_j) 的函数值之和,假设我们已经做完了第一部分得到了 (g(i)) 表示 (i) 以内质数的函数值之和,那么就这样递归:

[g(n)-sp(j)+sum_{i,e}f(p_i^e) imes (S(n/p_i^e,i)+[e ot=1]) ]

后面的 ([e ot=1]) 就表示如果不是只提了一次,那么就是某个质数大于 (1) 的幂次,这部分的答案是我们没有统计到的,所以加上。这一部分的复杂度我不太清楚,所以跟杜教筛一样玄学啊

例题

[模板]Min_25筛

这道题就是把原函数值拆成 (x^2,x^3) 的完全积性函数 (g),其他就是正常写法啦。

今天写 ( t Min25) 的时候遇到了问题,一开始筛 (g) 的时候要注意只能这么写 p[i]*p[i]<=w[j] ,必须保证有合数可以筛要不然会出问题

#include <cstdio>
#include <iostream>
#include <cmath>
using namespace std;
const int N = 200000;
const int M = N+5;
const int MOD = 1e9+7;
const int inv = 166666668;
#define int long long
int read()
{
	int x=0,f=1;char c;
	while((c=getchar())<'0' || c>'9') {if(c=='-') f=-1;}
	while(c>='0' && c<='9') {x=(x<<3)+(x<<1)+(c^48);c=getchar();}
	return x*f;
}
int n,cnt,tot,p[M],vis[M],w[M],sp1[M],sp2[M];
int sqr,id1[M],id2[M],g1[M],g2[M];
void init(int n)
{
	for(int i=2;i<=n;i++)
	{
		if(!vis[i])
		{
			p[++cnt]=i;
			sp1[cnt]=(sp1[cnt-1]+i)%MOD;
			sp2[cnt]=(sp2[cnt-1]+i*i)%MOD;
		}
		for(int j=1;j<=cnt && i*p[j]<=n;j++)
		{
			vis[i*p[j]]=1;
			if(i%p[j]==0) break;
		}
	}
}
int S(int x,int y)
{
	if(x<=p[y]) return 0;
	int k=x<=sqr?id1[x]:id2[n/x];
	int ans=(g2[k]-g1[k]-(sp2[y]-sp1[y]))%MOD;
	for(int i=y+1;i<=cnt && p[i]*p[i]<=x;i++)
	{
		int pr=p[i];
		for(int e=1;pr<=x;e++,pr*=p[i])
		{
			int xx=pr%MOD;
			ans=(ans+xx*(xx-1)%MOD*(S(x/pr,i)+(e!=1)))%MOD;
		}
	}
	return ans;
}
signed main()
{
	init(N);
	n=read();sqr=sqrt(n);
	for(int l=1,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		w[++tot]=n/l;
		g1[tot]=w[tot]%MOD;//要模要模 
		g2[tot]=g1[tot]*(g1[tot]+1)%MOD*(2*g1[tot]+1)%MOD*inv%MOD-1;
		g1[tot]=g1[tot]*(g1[tot]+1)/2%MOD-1;
		if(n/l<=sqr) id1[n/l]=tot;//编号方式是x 
		else id2[n/(n/l)]=tot;//编号方式是n/x
	}
	for(int i=1;i<=cnt;i++)
		for(int j=1;j<=tot && p[i]*p[i]<=w[j];j++)
		{
			int k=w[j]/p[i]<=sqr?id1[w[j]/p[i]]:id2[n/(w[j]/p[i])];
			g1[j]=(g1[j]-p[i]*(g1[k]-sp1[i-1]))%MOD;
			g2[j]=(g2[j]-p[i]*p[i]%MOD*(g2[k]-sp2[i-1]))%MOD;
		}
	printf("%lld
",(S(n,0)+1+MOD)%MOD);
}
原文地址:https://www.cnblogs.com/C202044zxy/p/14355392.html