博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
jzoj5683. 【GDSOI2018模拟4.22】Prime (Min_25筛+拉格朗日插值+主席树)
阅读量:6824 次
发布时间:2019-06-26

本文共 3203 字,大约阅读时间需要 10 分钟。

题面

1442599-20190221200232296-418723455.png

\(n\leq 10^{12},k\leq 100\)

题解

一眼就是一个\(Min\_25\)筛+拉格朗日插值优化,然而打完之后交上去发现只有\(60\)

\(tm\)还要用主席树优化……

大概是这样,设\(g(n,j)\)表示\(1\)\(n\)之间的所有满足\(i\)是质数或者\(i\)的最小质因子大于\(p_j\)的所有\(f(i)\)之和,我们根据递归地来求解\(g\),设一个阈值\(L=\sqrt{n}\),当\(n\leq L\)的时候,用主席树优化,能做到每一次查询只有\(O(\log\sqrt{n})\)

发现这玩意儿和杜教筛很像,所以就当它的复杂度是\(O(n^{\frac{2}{3}})\)好了……我实在算不来……

然而交上去还是\(T\)……你这才发现这题时间和空间都卡得要命……那么只好卡常了……

1.阈值\(L\)开到\(2100000\)左右实测最优,尽量不要开小,会需要多算很多东西,开大的话空间又会不够

2.空间卡着开,不需要多开的就别开了

3.在拉格朗日差值计算\(\sum_{i=2}^mf(i)\)的前缀和的时候,因为涉及到的所有\(m\)只有\(O(\sqrt{n})\)个,对于小于\(L\)的可以打表预处理,对于大于\(L\)\(m\)我们可以计算之后存起来。设\(x=\left\lfloor\frac{n}{m}\right\rfloor\),容易发现这里的\(x\leq \sqrt{n}\),想想整除分块,发现这里的\(x\)就是最大的满足\(\left\lfloor\frac{n}{x}\right\rfloor=m\)的数,那么每一个\(x\)都和\(m\)是一一对应的,我们就可以把它给存起来

4.还是过不去,吸氧好了

//minamoto#include
#define R register#define ll long long#pragma GCC optimize(3)#define fp(i,a,b) for(R int i=a,I=b+1;i
I;--i)#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)const int N=2.1e6+5,M=N*18,P=1e9+7;inline int add(R int x,R int y){return x+y>=P?x+y-P:x+y;}inline int dec(R int x,R int y){return x-y<0?x-y+P:x-y;}inline int mul(R int x,R int y){return 1ll*x*y-1ll*x*y/P*P;}int ksm(R int x,R int y){ R int res=1; for(;y;y>>=1,x=mul(x,x))if(y&1)res=mul(res,x); return res;}int p[N/3],sp[N/3],vis[N],f[N],sum[N],mnp[N],ssr[N],ss[N];int rt[N],ls[M],rs[M],s[M],fs[109],inv[109];int k,tot,m,sqr,cnt;ll n;void ins(int &p,int q,int l,int r,int x,int v){ s[p=++cnt]=add(s[q],v);if(l==r)return; int mid=(l+r)>>1; x<=mid?(rs[p]=rs[q],ins(ls[p],ls[q],l,mid,x,v)):(ls[p]=ls[q],ins(rs[p],rs[q],mid+1,r,x,v));}int query(int p,int l,int r,int x){ if(!p||l==r)return s[p]; int mid=(l+r)>>1; return x<=mid?add(s[rs[p]],query(ls[p],l,mid,x)):query(rs[p],mid+1,r,x);}void init(int n){ f[0]=sum[0]=P-1; fp(i,2,n){ if(!vis[i])p[++tot]=i,ssr[i]=f[i]=ksm(i,k),sp[tot]=add(sp[tot-1],f[i]); for(R int j=1;j<=tot&&1ll*i*p[j]<=n;++j){ vis[i*p[j]]=1,f[i*p[j]]=mul(f[i],f[p[j]]),mnp[i*p[j]]=j; if(i%p[j]==0)break; } } fp(i,2,n){ sum[i]=add(sum[i-1],f[i]),ssr[i]=add(ssr[i-1],ssr[i]); if(rt[i]=rt[i-1],mnp[i])ins(rt[i],rt[i-1],1,tot,mnp[i],f[i]); } inv[0]=inv[1]=1;fp(i,2,k+5)inv[i]=1ll*(P-P/i)*inv[P%i]%P;}int Lar(ll g,int n){ int k=g%P; if(k<=sqr)return sum[k]; if(ss[::n/g])return ss[::n/g]; int res=0,tmp=1,ty=(n&1)?P-1:1; fp(i,1,n)tmp=mul(tmp,inv[i]); fs[n+1]=1;fd(i,n,0)fs[i]=mul(fs[i+1],k-i); fp(i,0,n){ res=add(res,1ll*tmp*ty%P*sum[i]%P*fs[i+1]%P), tmp=1ll*tmp*(k-i)%P*(n-i)%P*inv[i+1]%P, ty=P-ty; }return ss[::n/g]=res;}int G(ll n,int m){ if(n<=sqr&&1ll*p[m+1]*p[m+1]>n)return ssr[n]; if(n<=sqr)return add(query(rt[n],1,tot,m+1),ssr[n]); if(!m)return Lar(n,k+1); while(1ll*p[m]*p[m]>n)--m; int res=Lar(n,k+1); while(m)res=dec(res,1ll*f[p[m]]*dec(G(n/p[m],m-1),sp[m-1])%P),--m; return res;}int main(){// freopen("testdata.in","r",stdin); freopen("prime.in","r",stdin); freopen("prime.out","w",stdout); scanf("%lld%d",&n,&k),sqr=2.1e6,init(sqr); printf("%d\n",G(n,tot-1)); return 0;}

转载于:https://www.cnblogs.com/bztMinamoto/p/10415037.html

你可能感兴趣的文章
struts2请求过程源码分析
查看>>
效率比较--集合
查看>>
jmeter IF控制器学习--使用实例
查看>>
memory prefix retro,re out 2
查看>>
WebDriver API 实例详解(四)
查看>>
dom01
查看>>
Android实例-如何使用系统剪切板(XE8+小米2)
查看>>
BAT-显示桌面图标
查看>>
PDO vs. MySQLi 选择哪一个?(PDO vs. MySQLi: Which Should You Use?)-转载
查看>>
信息安全系统设计基础第七周总结
查看>>
创建线程的三种方式
查看>>
Android项目依赖于第三方库(非jar包文件)
查看>>
cas HttpServletRequestWrapperFilter
查看>>
【Javascript第二重境界】函数
查看>>
SpringBoot 与 Web开发
查看>>
JavaWeb 三层框架
查看>>
BOOL, BOOLEAN, bool
查看>>
Mac 下 SVN 的使用
查看>>
简述session
查看>>
Android APK反编译教程(带工具)
查看>>