标签：prime 25 min primeorminp 质数 minp 详解 pj quad

扯淡

min_25筛是由min_25提出的求积性函数前缀和的亚线性算法，和一个叫“扩展埃氏筛”的东西有着微妙的关系。

至于是什么关系，我也不太清楚，反正有人说很像有人说就是一个东西（雾）

这段话并不是废话

约定

为了方便后面描述，这里写一些用到的约定和符号表示，以免产生恐惧

$1$1被开除正整数籍 也就是说“前缀和”之类的都是从$2$2开始，对答案所求的前缀和同样，最后手动加$1$1

$\pi(n)$π(n)表示$1 \sim n$1∼n中质数个数的规模（其实问题不大，后面就懂了）

$p_i$pi​表示第$i$i个质数，单独的$p$p均表示质数

$prime$prime表示质数集合 $minp(i)$minp(i)表示$i$i的最小质因子

$p^c$pc表示一个只含一个质因子的数

流程

首先所求的函数$f(n)$f(n)需要满足：

1. 是个积性函数
2. 在质数处的取值$f(p)$f(p)是一个关于$p$p的多项式。我们把这个多项式拆成若干单项式分别计算在相加，这样就变成了"$f(p)$f(p)的值是一个关于$p$p的单项式"，我们记为$f(p)=p^k$f(p)=pk
3. $f(p^c)$f(pc)的值可以快速计算

首先考虑计算这个东西

$\sum_{i=2}^n[i\in prime]i^k$i=2∑n​[i∈prime]ik

注意是$i^k$ik不是$f(i)$f(i)，虽然这里还没有区别，但它们只是质数位置相等

首先既然叫扩展埃氏筛先回忆一下埃氏筛在干什么：开始写下所有正整数，然后从小到大，如果一个没被筛说明它是质数，用它筛掉后面的倍数

我们可以用这个思路计算上面的式子，先算出所有的和，再用质数筛掉所有合数项

设$g(n,j)$g(n,j)表示用前$j$j个质数筛了之后剩余项的和

$g(n,j)=\sum_{i=2}^n[i \in prime \quad or\quad minp(i)>p_j]i^k$g(n,j)=i=2∑n​[i∈primeorminp(i)>pj​]ik

和就是$g(n,\pi(\sqrt{n}))$g(n,π(n​))

注意合数的最小质因子不会超过$\sqrt n$n​,所以如果$p_j^2>n$pj2​>n有$g(n,j)=g(n,j-1)$g(n,j)=g(n,j−1)

对剩下情况考虑转移，即从$g(n,j-1)$g(n,j−1)去掉被$p_j$pj​删掉的项

$g(n,j-1)=\sum_{i=2}^n[i\in prime \quad or \quad minp(i)\geq p_j]$g(n,j−1)=i=2∑n​[i∈primeorminp(i)≥pj​]

如果脑子转不过来，可以强行算

$[i\in prime \quad or \quad minp(i)\geq p_j]$[i∈primeorminp(i)≥pj​]且不满足$[i \in prime \quad or\quad minp(i)>p_j]$[i∈primeorminp(i)>pj​]
即
$[i\in prime \quad or \quad minp(i)\geq p_j]$[i∈primeorminp(i)≥pj​]且$[i \notin prime] 且 [minp(i)\leq p_j]$[i∈/​prime]且[minp(i)≤pj​]
综上
$[i \notin prime,minp(i)=p_j]$[i∈/​prime,minp(i)=pj​]

所以

$g(n,j)=g(n,j-1)-\sum_{i=2}^n[i \notin prime,minp(i)=p_j]i^k$g(n,j)=g(n,j−1)−i=2∑n​[i∈/​prime,minp(i)=pj​]ik

可以提一个$p_j$pj​出来,因为不能有质数，就把$p_j$pj​去掉，刚好是$1$1

$g(n,j)=g(n,j-1)-p_j^k\sum_{i=2}^{\lfloor\frac{n}{p_j}\rfloor}[minp(i)\geq p_j]i^k$g(n,j)=g(n,j−1)−pjk​i=2∑⌊pj​n​⌋​[minp(i)≥pj​]ik

再次观察

$g(n,j)=\sum_{i=2}^n[i \in prime \quad or\quad minp(i)>p_j]i^k$g(n,j)=i=2∑n​[i∈primeorminp(i)>pj​]ik

发现可以拆成小于$p_j$pj​的质数和大于等于$p_j$pj​的所有数，第二个就是刚才的式子

$g(n,j)=g(n,j-1)-p_j^k(g(\lfloor\frac{n}{p_j}\rfloor,j-1)-\sum_{i=1}^{j-1}p_i^k)$g(n,j)=g(n,j−1)−pjk​(g(⌊pj​n​⌋,j−1)−i=1∑j−1​pik​)

注意$n$n由于都是一直整除，所以只会有$O(\sqrt n)$O(n​)种取值，可以整除分块找出来强行离散化。在记录某个值$v$v所在位置的时候，如果$v>\sqrt n$v>n​,我们另开一个数组存到$\lfloor\frac{n}{v}\rfloor$⌊vn​⌋里面

然后后面只会用到最后一项，所以第二维可以滚掉

求答案时，设

$S(n,j)=\sum_{i=2}^n[minp(i)>p_j]f(i)$S(n,j)=i=2∑n​[minp(i)>pj​]f(i)

分质数和合数分别计算

质数部分用$g$g算出来的减去前$j$j个

$g(n,\pi(\sqrt n))-\sum_{i=1}^jp_i$g(n,π(n​))−i=1∑j​pi​

合数部分枚举最小质因子和它的次数

$\sum_{k=j+1}^{p_k\leq\sqrt n}\sum_{e=1}^{p_k^e\leq n}f(p_k^e)(S(\lfloor\frac{n}{p_k^e}\rfloor,k)+[e>1])$k=j+1∑pk​≤n​​e=1∑pke​≤n​f(pke​)(S(⌊pke​n​⌋,k)+[e>1])

$[e>1]$[e>1]指如果指数大于$1$1它本身就是合数，需要统计答案

两部分相加即可

总复杂度大概$O(n^{2\over3})$O(n32​)常数较小，大约$2s$2s过$1e10$1e10,$3s$3s过$1e11$1e11

模板题

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#include <cmath>

#define MAXN 200005
using namespace std;
typedef long long ll;
const int MOD=1e9+7,INV6=(MOD+1)/6;
inline int add(const int& x,const int& y){return x+y>=MOD? x+y-MOD:x+y;}
inline int dec(const int& x,const int& y){return x<y? x-y+MOD:x-y;}
int np[MAXN],pl[MAXN],cnt;
inline void init(const int& N)
{
	np[1]=1;
	for (int i=2;i<=N;i++)
	{
		if (!np[i]) pl[++cnt]=i;
		int x;
		for (int j=1;(x=i*pl[j])<=N;j++)
		{
			np[x]=1;
			if (i%pl[j]==0) break;
		}
	}
}
ll val[MAXN];
int tot;
int g1[MAXN],g2[MAXN],sum1[MAXN],sum2[MAXN];
int key[MAXN],yek[MAXN];
ll n;
int m;
inline int getkey(const ll& v){return v<=m? key[v]:yek[n/v];}
int S(ll n,int j)
{
	if (pl[j]>=n) return 0;
	int k=getkey(n);
	int ans=dec(dec(g2[k],g1[k]),dec(sum2[j],sum1[j]));
	for (int k=j+1;k<=cnt&&(ll)pl[k]*pl[k]<=n;k++)
		for (ll e=1,pe=pl[k];pe<=n;e++,pe*=pl[k])
			ans=add(ans,pe%MOD*(pe%MOD-1)%MOD*(S(n/pe,k)+(e>1))%MOD);				
	return ans;
}
int main()
{
	scanf("%lld",&n);
	m=sqrt(n);
	init(m);
	for (ll l=1,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		val[++tot]=n/l;
		int t=val[tot]%MOD;
		g1[tot]=((ll)t*(t+1)/2)%MOD-1;
		g2[tot]=(ll)t*(t+1)%MOD*(2*t+1)%MOD*INV6%MOD-1;
		if (val[tot]<=m) key[val[tot]]=tot;
		else yek[n/(n/l)]=tot;
	}
	for (int j=1;j<=cnt;j++)
	{
		for (int i=1;i<=tot&&(ll)pl[j]*pl[j]<=val[i];i++)
		{
			int k=getkey(val[i]/pl[j]);
			g1[i]=dec(g1[i],(ll)pl[j]*dec(g1[k],sum1[j-1])%MOD);
			g2[i]=dec(g2[i],(ll)pl[j]*pl[j]%MOD*dec(g2[k],sum2[j-1])%MOD);
		}		
		sum1[j]=add(sum1[j-1],pl[j]),sum2[j]=add(sum2[j-1],(ll)pl[j]*pl[j]%MOD);
	}
	printf("%d\n",S(n,0)+1);
	return 0;
}

• 点赞
• 收藏
• 分享
• • 文章举报

Lstdo 发布了56 篇原创文章 · 获赞 3 · 访问量 4765 私信 关注

标签：prime,25,min,primeorminp,质数,minp,详解,pj,quad
来源： https://blog.csdn.net/luositing/article/details/104005161

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