LGP7884题解

是的,这是一篇使用 min25 筛的题解。。。

本题解参考command_block大佬的博客,代码是对其在 LOJ 上的提交卡常后写出来的。

ML 板子把数据开到 \(10^{13}\) 速度还和供题人的 ML 速度差不多快就离谱。。。好吧我吸了氧

这个板子的原理是使用树状数组优化的 min25 筛,在下面会详细讲解。复杂度是 \(O((\frac n {\log n})^{\frac 2 3})\) 的。

首先有经典 DP:

\[G(n,k)=G(n,k-1)-G(\frac n {pri_k},k-1)+G(pri_{k-1},k-1)) \]

边界条件为:

\[G(n,0)=n \]

根据积分,可以得到朴素的 DP 是 \(O(\frac {n^{\frac 3 4}} {\log n})\) 的。接下来尝试使用树状数组维护对于 \(n\) 过小时的 \(G\)

设分治线为 \(B1\)。对 \(B1\) 一下的所有 \(G\) 使用树状数组维护,对于 \(B1\) 以上的 \(G\) 使用 DP。

这一部分的复杂度为 \(O(B1\log B1+\frac n {\sqrt{B1}\log n})\),取 \(B1=\frac {n^{\frac 2 3}} {\log^{\frac 4 3} n}\) 可以得到复杂度为 \(O(\frac {n^{\frac 2 3}} {\log^{\frac 1 3} n})\)

继续分治,令 \(B2\) 以下的部分暴力统计。

\(B2=\sqrt [6] n\) 时,复杂度为 \(O(n^{\frac 2 3}\log n+\frac n {\sqrt{B1}\log n}+B1)\),取 \(B1=(\frac n {\log n})^{\frac 2 3}\) 可以得到复杂度为 \(O((\frac n {\log n})^{\frac 2 3})\)

复杂度的具体推导可以看 blog。笔者看不懂于是懒得解释直接贺代码了

//感谢@command_block的板子 
#include<cstdio>
#include<cmath>
const int M=3e7,Lim=8.5e7+10;
typedef unsigned long long ull;
int l2,tot,lim,BIT[Lim];
ull n,g[M],w[M];
double inv[M];
bool e[Lim];
inline ull min(const ull&a,const ull&b){
	return a>b?b:a;
}
inline ull max(const ull&a,const ull&b){
	return a>b?a:b;
}
inline void Add(int x){
	e[x]=1;
	while(x<=lim)++BIT[x],x+=x&-x;
}
inline ull Query(int x){
	ull sum=x;
	while(x)sum-=BIT[x],x^=x&-x;
	return sum;
}
signed main(){
	register int i,j,tl,tl2,tl3;
	register ull t,r,x0;
	scanf("%llu",&n);
	lim=min(max(max((pow(n/log2(n),0.66)),l2=sqrtl(n)+1),10000),n);
	for(i=1;i<=l2;i++)w[i]=i-1,inv[i]=1./i;
	for(tot=1;1ull*lim*tot<n;tot++)g[tot]=n*inv[tot]+1e-6-1;
	--tot;Add(1);
	for(i=2;1ull*i*i<=n;i++){
		if(e[i])continue;
		x0=w[i-1];t=n/i;r=1ull*i*i;
		tl=min(n/r,(ull)tot);tl2=min(tl,n/(1ull*l2*i));tl3=min(tl2,tot/i);
		for(j=1;j<=tl3;++j)g[j]-=g[j*i]-x0;
		for(j=tl3+1;j<=tl2;j++)g[j]-=Query(t*inv[j]+1e-6)-x0;
		for(j=tl2+1;j<=tl;++j)g[j]-=w[int(t*inv[j]+1e-6)]-x0;
		for(j=l2;j>=r;--j)w[j]-=w[int(j*inv[i]+1e-6)]-x0;
		if(1ull*i*i<=lim){
			for(j=i*i;j<=lim;j+=i)if(!e[j])Add(j);
		}
	}
	if(!tot)g[1]=Query(n);
	printf("%llu",g[1]);
}
posted @ 2022-01-11 14:05  Prean  阅读(18)  评论(0编辑  收藏  举报
var canShowAdsense=function(){return !!0};