【XSY3473】Frog(min25筛)

一些记号:

  • P \mathbb{P} P 为质数集, p i p_i pi 表示第 i i i 个质数。
  • lpf ⁡ ( x ) \operatorname{lpf}(x) lpf(x) 表示 x x x 的最小质因数为第几个质数。
  • 特别地,令 lpf ⁡ ( 1 ) = ∞ \operatorname{lpf}(1)=\infty lpf(1)=

c i c_i ci 表示 p i p_i pi m ! m! m! 中的次数,注意到当 p i > m p_i>\sqrt m pi>m 时, c i = ⌊ m p i ⌋ c_i=\lfloor\dfrac{m}{p_i}\rfloor ci=pim,所以 c i c_i ci 的取值只有 O ( m ) O(\sqrt m) O(m ) 种,而且 c i c_i ci 单调不降。

F ( i , n ) F(i,n) F(i,n) 表示不考虑前 i i i 个质数时的答案,即:
F ( i , n ) = ∑ x = 1 n ∑ d ∣ m ! [ lpf ⁡ ( x d ) > i ] σ 0 ( x d ) F(i,n)=\sum_{x=1}^n \sum_{d|m!}[\operatorname{lpf}(xd)>i]\sigma_0(xd) F(i,n)=x=1ndm![lpf(xd)>i]σ0(xd)
注意 x d = 1 xd=1 xd=1 的情况也会被算入。

考虑递推到 F ( i − 1 , n ) F(i-1,n) F(i1,n),考虑 p i p_i pi 分别在 x x x d d d 中出现的次数:
F ( i − 1 , n ) = ∑ k ≥ 0 , p i k ≤ n ∑ 0 ≤ j ≤ c i ( k + j + 1 ) F ( i , ⌊ n p i k ⌋ ) = ∑ k ≥ 0 , p i k ≤ n 1 2 ( c i + 1 ) ( c i + 2 k + 2 ) F ( i , ⌊ n p i k ⌋ ) \begin{aligned} F(i-1,n)&=\sum_{k\geq 0,p_i^k\leq n}\sum_{0\leq j\leq c_i}(k+j+1) F(i,\lfloor\frac{n}{p_i^k}\rfloor)\\ &=\sum_{k\geq 0,p_i^k\leq n}\frac{1}{2}(c_i+1)(c_i+2k+2)F(i,\lfloor\frac{n}{p_i^k}\rfloor) \end{aligned} F(i1,n)=k0,pikn0jci(k+j+1)F(i,pikn)=k0,pikn21(ci+1)(ci+2k+2)F(i,pikn)
如果我们能处理出 p i + 1 2 > n p_{i+1}^2>n pi+12>n 时的边界情况,这个递推的时间复杂度就和 min25 筛是一样的了。

p i + 1 2 > n p_{i+1}^2>n pi+12>n 时,直接考虑定义式, F ( i , n ) F(i,n) F(i,n) 即为:
F ( i , n ) = ∑ d ∣ m ! [ lpf ⁡ ( d ) > i ] σ 0 ( d ) + ∑ j > i , p j ≤ n ∑ d ∣ m ! [ lpf ⁡ ( d ) > i ] σ 0 ( p j d ) \begin{aligned} F(i,n)&=\sum_{d|m!}[\operatorname{lpf}(d)>i]\sigma_0(d)+\sum_{j>i,p_j\leq n}\sum_{d|m!}[\operatorname{lpf}(d)>i]\sigma_0(p_jd) \end{aligned} F(i,n)=dm![lpf(d)>i]σ0(d)+j>i,pjndm![lpf(d)>i]σ0(pjd)
我们先考虑一个子问题 ∑ d ∣ T σ 0 ( d ) \sum\limits_{d|T}\sigma_0(d) dTσ0(d),其中设 T = ∏ i p i t i T=\prod\limits_i p_i^{t_i} T=ipiti

相当于对每一个质因数 p i p_i pi 选一个 t i ′ ≤ t i t'_i\leq t_i titi 以确定 d d d,然后这个 d d d 的贡献是 ∏ i ( 1 + t i ′ ) \prod\limits_i (1+t_i') i(1+ti)。于是:
∑ d ∣ T σ 0 ( d ) = ∏ i ( ∑ t i ′ = 0 t i ( 1 + t i ′ ) ) = ∏ i 1 2 ( t i + 1 ) ( t i + 2 ) \begin{aligned} \sum\limits_{d|T}\sigma_0(d)&=\prod_{i}\left(\sum_{t_i'=0}^{t_i}(1+t_i')\right)\\ &=\prod_i \frac{1}{2}(t_i+1)(t_i+2) \end{aligned} dTσ0(d)=iti=0ti(1+ti)=i21(ti+1)(ti+2)
然后再考虑 ∑ d ∣ T σ 0 ( p j d ) \sum\limits_{d|T}\sigma_0(p_jd) dTσ0(pjd)
∑ d ∣ T σ 0 ( p j d ) = ∑ d ∣ ( p j T ) σ 0 ( d ) − ∑ d ∣ ( T / p j t j ) σ 0 ( d ) = ∑ d ∣ T σ 0 ( d ) 1 2 ( t j + 2 ) ( t j + 3 ) − 1 1 2 ( t j + 1 ) ( t j + 2 ) \begin{aligned} \sum_{d|T}\sigma_0(p_jd)&=\sum_{d|(p_jT)}\sigma_0(d)-\sum_{d|(T/p_j^{t_j})}\sigma_0(d)\\ &=\sum_{d|T}\sigma_0(d)\frac{\frac{1}{2}(t_j+2)(t_j+3)-1}{\frac{1}{2}(t_j+1)(t_j+2)}\\ \end{aligned} dTσ0(pjd)=d(pjT)σ0(d)d(T/pjtj)σ0(d)=dTσ0(d)21(tj+1)(tj+2)21(tj+2)(tj+3)1
再看回原式, ∑ d ∣ m ! [ lpf ⁡ ( d ) > i ] σ 0 ( p j d ) \sum\limits_{d|m!}[\operatorname{lpf}(d)>i]\sigma_0(p_jd) dm![lpf(d)>i]σ0(pjd) 其实可以转化成 ∑ d ∣ ( m ! / ∏ k ≤ i p k c k ) σ 0 ( p j d ) \sum\limits_{d|(m!/\prod_{k\leq i}p_k^{c_k})}\sigma_0(p_jd) d(m!/kipkck)σ0(pjd),于是有:
F ( i , n ) = ∑ d ∣ m ! [ lpf ⁡ ( d ) > i ] σ 0 ( d ) + ∑ j > i , p j ≤ n ∑ d ∣ m ! [ lpf ⁡ ( d ) > i ] σ 0 ( p j d ) = ∑ d ∣ ( m ! / ∏ k ≤ i p k c k ) σ 0 ( d ) + ∑ j > i , p j ≤ n ∑ d ∣ ( m ! / ∏ k ≤ i p k c k ) σ 0 ( p j d ) = ∑ d ∣ ( m ! / ∏ k ≤ i p k c k ) σ 0 ( d ) + ∑ j > i , p j ≤ n ∑ d ∣ ( m ! / ∏ k ≤ i p k c k ) σ 0 ( d ) 1 2 ( c j + 2 ) ( c j + 3 ) − 1 1 2 ( c j + 1 ) ( c j + 2 ) = ( ∏ j > i 1 2 ( c i + 1 ) ( c i + 2 ) ) ( 1 + ∑ j > i , p j ≤ n 1 2 ( c j + 2 ) ( c j + 3 ) − 1 1 2 ( c j + 1 ) ( c j + 2 ) ) \begin{aligned} F(i,n)&=\sum_{d|m!}[\operatorname{lpf}(d)>i]\sigma_0(d)+\sum_{j>i,p_j\leq n}\sum_{d|m!}[\operatorname{lpf}(d)>i]\sigma_0(p_jd)\\ &=\sum\limits_{d|(m!/\prod_{k\leq i}p_k^{c_k})}\sigma_0(d)+\sum_{j>i,p_j\leq n}\sum\limits_{d|(m!/\prod_{k\leq i}p_k^{c_k})}\sigma_0(p_jd)\\ &=\sum\limits_{d|(m!/\prod_{k\leq i}p_k^{c_k})}\sigma_0(d)+\sum_{j>i,p_j\leq n}\sum_{d|(m!/\prod_{k\leq i}p_k^{c_k})}\sigma_0(d)\frac{\frac{1}{2}(c_j+2)(c_j+3)-1}{\frac{1}{2}(c_j+1)(c_j+2)}\\ &=\left(\prod_{j>i}\frac{1}{2}(c_i+1)(c_i+2)\right)\left(1+\sum_{j>i,p_j\leq n}\frac{\frac{1}{2}(c_j+2)(c_j+3)-1}{\frac{1}{2}(c_j+1)(c_j+2)}\right) \end{aligned} F(i,n)=dm![lpf(d)>i]σ0(d)+j>i,pjndm![lpf(d)>i]σ0(pjd)=d(m!/kipkck)σ0(d)+j>i,pjnd(m!/kipkck)σ0(pjd)=d(m!/kipkck)σ0(d)+j>i,pjnd(m!/kipkck)σ0(d)21(cj+1)(cj+2)21(cj+2)(cj+3)1=(j>i21(ci+1)(ci+2))1+j>i,pjn21(cj+1)(cj+2)21(cj+2)(cj+3)1
相当于我在数轴上有 O ( m ) O(\sqrt m) O(m ) 个关键点,关键点之间的质数的 c c c 都相等,我在数轴上还有 O ( n ) O(\sqrt n) O(n ) 个关键点用于询问,于是我们把这 O ( n ) + O ( m ) O(\sqrt n)+O(\sqrt m) O(n )+O(m ) 个关键点混在一起用 min25 筛出来它们的素数个数即可。

有关 min25 筛复杂度的证明:

考虑 F ( i , j ) F(i,j) F(i,j) 的状态数,第二维的取值为 n , ⌊ n 2 ⌋ , ⋯   , ⌊ n n ⌋ , n , ⋯   , 1 n,\lfloor\frac{n}{2}\rfloor,\cdots,\lfloor\frac{n}{\sqrt n}\rfloor,\sqrt n,\cdots,1 n,2n,,n n,n ,,1 O ( n ) O(\sqrt n) O(n ) 种(整除分块),第一维的取值范围是 O ( π ( j ) ) O(\pi(\sqrt j)) O(π(j )),于是共有状态数:
∑ i = 1 n O ( π ( i ) ) + ∑ i = 1 n O ( π ( n i ) ) ≈ ∑ i = 1 n O  ⁣ ( i ln ⁡ i ) + ∑ i = 1 n O  ⁣ ( n i ln ⁡ n i ) ≈ O ( ∫ 1 n x ln ⁡ x  d x + ∫ 1 n n x ln ⁡ n x  d x ) ≈ O ( n 3 4 ln ⁡ n ) \begin{aligned} &\sum_{i=1}^{\sqrt n}O(\pi(\sqrt i))+\sum_{i=1}^{\sqrt n}O(\pi(\sqrt{\frac{n}{i}}))\\ \approx&\sum_{i=1}^{\sqrt n}O\!\left(\frac{\sqrt i}{\ln \sqrt i}\right)+\sum_{i=1}^{\sqrt n}O\!\left(\frac{\sqrt{\frac{n}{i}}}{\ln \sqrt{\frac{n}{i}}}\right)\\ \approx&O\left(\int_1^{\sqrt n}\frac{\sqrt x}{\ln \sqrt x}\text{ d}x+\int_{1}^{\sqrt n}\frac{\sqrt{\frac{n}{x}}}{\ln \sqrt{\frac{n}{x}}}\text{ d}x\right)\\ \approx&O\left(\frac{n^{\frac{3}{4}}}{\ln n}\right) \end{aligned} i=1n O(π(i ))+i=1n O(π(in ))i=1n O(lni i )+i=1n O(lnin in )O(1n lnx x  dx+1n lnxn xn  dx)O(lnnn43)
转移时的复杂度忽略,否则太难分析了(

所以你可以把复杂度当成 O ( n 3 4 ) O(n^{\frac{3}{4}}) O(n43)

写代码时写 id 那一段时竟然绕着绕着给自己绕晕了,后来发现只需要记住一件事:id 只是一种给关键点的编码方式,所以你只需要找到一种编码方式保证不会存在两个关键点编号相同即可。

数论题代码果然及其难写,所以代码也写得及其丑,跑得也及其慢。看啥时候有空规范一下自己数论题的码风(

#include<bits/stdc++.h>

#define N 100010
#define ll long long
#define LNF 0x7fffffffffffffff

using namespace std;

namespace modular
{
	const int mod=323232323,inv2=(mod+1)/2;
	inline int add(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;}
	inline int mul(int x,int y){return 1ll*x*y%mod;}
	inline void Add(int &x,int y){x=x+y>=mod?x+y-mod:x+y;}
	inline void Dec(int &x,int y){x=x-y<0?x-y+mod:x-y;}
	inline void Mul(int &x,int y){x=1ll*x*y%mod;}
}using namespace modular;

inline int poww(int a,int b)
{
	int ans=1;
	while(b)
	{
		if(b&1) ans=mul(ans,a);
		a=mul(a,a);
		b>>=1;
	}
	return ans;
}

inline ll read()
{
	ll x=0;
	int f=1;
	char ch=getchar();
	while(ch<'0'||ch>'9')
	{
		if(ch=='-') f=-1;
		ch=getchar();
	}
	while(ch>='0'&&ch<='9')
	{
		x=(x<<1)+(x<<3)+(ch^'0');
		ch=getchar();
	}
	return x*f;
}

int tot;
int cnt,prime[N];
bool notprime[N];
ll n,m,sn,sm;
ll kp[N*5];

namespace ID
{
	int idn1[N],idn2[N],idm1[N],idm2[N],idmm[N];
	int &id(ll x)
	{
		if(x<=n&&x==n/(n/x)) return x<=sn?idn1[x]:idn2[n/x];
		if(x<=m&&x==m/(m/x)) return x<=sm?idm1[x]:idm2[m/x];
		return idmm[x];
	}
}using ID::id;

int is1[N],is2[N];

ll get(ll p)
{
	ll x=m,ans=0;
	while(x)
	{
		x/=p;
		ans+=x;
	}
	return ans;
}

ll c[N];

void init()
{
	const int nn=100000;
	is1[0]=1;
	for(int i=2;i<=nn;i++)
	{
		if(!notprime[i])
		{
			prime[++cnt]=i;
			c[cnt]=get(i);
			is1[cnt]=mul(is1[cnt-1],poww(mul(inv2,mul((c[cnt]+1)%mod,(c[cnt]+2)%mod)),mod-2));
			int a=dec(mul(inv2,mul((c[cnt]+2)%mod,(c[cnt]+3)%mod)),1);
			int b=mul(inv2,mul((c[cnt]+1)%mod,(c[cnt]+2)%mod));
			is2[cnt]=add(is2[cnt-1],mul(a,poww(b,mod-2)));
		}
		for(int j=1;j<=cnt&&i*prime[j]<=nn;j++)
		{
			notprime[i*prime[j]]=1;
			if(!(i%prime[j])) break;
		}
	}
}

ll g[N*5];

void work1()
{
	for(int i=1;i<=tot;i++) g[i]=kp[i]-1;
	for(int i=1;i<=cnt;i++)
	{
		ll p=prime[i],p2=p*p;
		for(int j=tot;p2<=kp[j];j--)
			g[j]-=g[id(kp[j]/p)]-(i-1);
	}
}

int f[N*5];
int sum1[N*5],sum2[N*5];
bool vis[N*5];

int bound(int i,ll n)
{
	int a=mul(sum1[tot],is1[i]);
	int b=prime[i]<=n?dec(sum2[id(n)],is2[i]):0;
	return mul(a,add(b,1));
}

void work2()
{
	for(int i=cnt;i>=1;i--)
	{
		ll p=prime[i],p2=p*p,p22=(i==cnt?LNF:1ll*prime[i+1]*prime[i+1]);
		for(int j=tot;p2<=kp[j];j--)
		{
			if(!vis[j]) f[j]=bound(i,kp[j]),vis[j]=1;
			Mul(f[j],mul(inv2,mul((c[i]+1)%mod,(c[i]+2)%mod)));
			ll pk=p;
			for(int k=1;pk<=kp[j];k++,pk*=p)
			{
				int coef=mul(inv2,mul((c[i]+1)%mod,(c[i]+2*k+2)%mod));
				ll v=kp[j]/pk;
				if(p22>v) Add(f[j],mul(coef,bound(i,v)));
				else Add(f[j],mul(coef,f[id(v)]));
			}
		}
	}
}

int main()
{
	n=read(),m=read();
	sn=sqrt(n),sm=sqrt(m);
	init();
	for(ll l=1,r;l<=n;l=r+1)
		kp[++tot]=r=(n/(n/l));
	for(ll l=1,r;l<=m;l=r+1)
		kp[++tot]=r=(m/(m/l));
	for(ll i=1;i*i<=m;i++)
		kp[++tot]=i;
	sort(kp+1,kp+tot+1);
	tot=unique(kp+1,kp+tot+1)-kp-1;
	for(int i=1;i<=tot;i++) id(kp[i])=i;
	work1();
	sum1[1]=1;
	//kp[i-1]+1 ~ kp[i]
	for(int i=2;i<=tot;i++)
	{
		ll np=g[i]-g[i-1];
		ll c=get(kp[i]);
		sum1[i]=mul(sum1[i-1],poww(mul(inv2,mul((c+1)%mod,(c+2)%mod)),np));
		int a=dec(mul(inv2,mul((c+2)%mod,(c+3)%mod)),1);
		int b=mul(inv2,mul((c+1)%mod,(c+2)%mod));
		sum2[i]=add(sum2[i-1],mul(np,mul(a,poww(b,mod-2))));
	}
	work2();
	printf("%d\n",f[id(n)]);
	return 0;
}
/*
10 3
*/
/*
323232323 23
*/
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值