ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

数论 筛法+拉格朗日插值 lgP5493题解

2022-01-10 16:32:31  阅读:152  来源: 互联网

标签:筛法 int 题解 ll register lgP5493 ans id mod


卡完常后来造福一下人类

如何从4.80s卡到920ms.jpg

本题解的复杂度为 \(O(\frac {n^{3/4}} {\log n})\),然而标算是 \(O(\frac {n^{2/3}} {\log^{1/3} n})\) 的。。。

有时间尝试卡一下标算,但是看样子好像已经卡过一些了,不知道能不能比我这个代码快(

首先亮出经典 DP:

\[f(n,id)=f(n,id-1)-p_{id}^k \times (f(\frac n {p_{id}},id-1)-f(p_{id},id-1)) \]

然后你写完之后稍微卡一下,再吸个氧就能得到4.80s的代码了。

稍微卡一下指把DP的部分中的int和ll分开,并且线性筛进行了一些神奇的优化(

然后我们开始卡。

首先加了一个FastMod,速度变成了2.57s

然后众所周知的是,实数除法比整数除法要快,变成了1.31s。

然后我们知道线性筛的原理是 用自身最小的质因子筛掉自己,那么我们没有必要用除法,将其记录下来即可,1.13s。

然后我们将减法优化改成暴力取模,发现变成了1.06s。

然后由于我的DP过程中边界是这样判的:

for(;j<=tot&&pri[i]<=(m1=w[j]*invp[i]);++j)

我们发现只需要将pri[top+1]改为INF就能够避免掉前面的那个j<=tot,这次卡进了1s,970ms。

然后由于我们DP时每次都计算了一遍sum[i-1]+p,我们就新开了一个变量s将其存下来,920ms。

upd:把f中的n故技重施能卡到915ms。

不知道还能不能卡/youl

upd:把前面的一些“优化”删掉之后跑了885ms。

#include<cstdio>
#include<cmath>
typedef long long ll;
typedef __uint128_t L;
typedef unsigned long long ull;
const ll M=2e5+5;
int k,p,a[15],ifac[15],sl[15],sr[15];
int S,id1[M],id2[M];ll tot,g[M<<1];ll n,w[M<<1];
int top,F[17985],pri[17985],sum[17985],pos[M];bool zhi[M];
double invp[17985];
struct FastMod{
    ull b,m;
    FastMod(ull b):b(b),m(ull((L(1)<<64)/b)){}
    friend inline ull operator%(const ull&a,const FastMod&mod){
        ull q=(L(mod.m)*a)>>64;
        ull r=a-q*mod.b;
        return r>=mod.b?r-mod.b:r;
    }
}mod(2);
inline int pow(int a,int b){
	register int ans=1;
	for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)ans=1ll*ans*a%mod;
	return ans;
}
inline void init(){
	register int i;k+=2;a[1]=sl[0]=sr[k+1]=ifac[0]=ifac[1]=1;
	for(i=2;i<=k;++i)a[i]=(a[i-1]+pow(i,k-2))%mod,ifac[i]=1ll*(p-p/i)*ifac[p%i]%mod;
	for(i=2;i<=k;++i)ifac[i]=1ll*ifac[i-1]*ifac[i]%mod;
	for(i=1;i<=k;++i)a[i]=1ll*ifac[i-1]*(k-i&1?p-ifac[k-i]:ifac[k-i])%mod*a[i]%mod;
}
inline int f(const int&n){
	register int i,N=n+p;register ull ans=0;
	for(i=1;i<=k;++i)sl[i]=1ll*sl[i-1]*(N-i)%mod;
	for(i=k;i>=1;--i)sr[i]=1ll*sr[i+1]*(N-i)%mod;
	for(i=1;i<=k;++i)ans+=1ll*sr[i+1]*sl[i-1]%mod*a[i];
	return ans%mod;
}
inline void sieve(const int&n){
	register int i=6,j,x,m;top=2;
	F[1]=pow(pri[1]=2,k);F[2]=pow(pri[2]=3,k);
	sum[1]=F[1];sum[2]=F[1]+F[2];
	invp[1]=1./2*(1+1e-15);invp[2]=1./3*(1+1e-15);
	do{
		if(!zhi[m=i-1]&&i-1<=n){
			pri[++top]=m;sum[top]=(sum[top-1]+(F[top]=pow(m,k)))%mod;
			invp[top]=1./m*(1+1e-15);
		}
		for(j=3;j<=top&&(x=m*pri[j])<=n;++j){
			zhi[x]=true;if((pos[x]=j)==pos[m])break;
		}
		if(!zhi[m=i+1]&&i+1<=n){
			pri[++top]=m;sum[top]=(sum[top-1]+(F[top]=pow(m,k)))%mod;
			invp[top]=1./m*(1+1e-15);
		}
		for(j=3;j<=top&&(x=m*pri[j])<=n;++j){
			zhi[x]=true;if((pos[x]=j)==pos[m])break;
		}
	}while((i+=6)-1<=n);pri[++top]=p;invp[top]=0;
}
void Solve(const ll&n){
	const ll&n9=n/1e9;
	register int i,j,k,s;register ll m,L=1,R;
	for(;L<=n;L=R+1,--g[tot]){
		R=n/(m=w[++tot]=1.*n/L);g[(m<=S?id1[m]:id2[R])=tot]=f(m%mod);
	}
	for(i=1;i<=top;++i){
		s=sum[i-1]+p;
		for(j=1;pri[i]<=(m=w[j]*invp[i]);++j){
			g[j]+=1ll*F[i]*(s-g[m<=S?id1[m]:id2[int(1.*n/m)]])%mod;
			if(g[j]>=p)g[j]-=p;
		}
	}
}
signed main(){
	register int i=1;register ull ans=0;
	scanf("%lld%d%d",&n,&k,&p);mod=FastMod(p);
	sieve(S=sqrt(n));init();Solve(n);
	for(register ll m;i<=S;++i){
		m=1.*n/i;
		ans+=1ll*i*i%mod*g[m<=S?id1[m]:id2[ll(1.*n/m)]]%mod;
		if(ans>=p)ans-=p;
	}
	printf("%d",ans);
}

标签:筛法,int,题解,ll,register,lgP5493,ans,id,mod
来源: https://www.cnblogs.com/lmpp/p/15784940.html

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

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有