• jzoj5683. 【GDSOI2018模拟4.22】Prime (Min_25筛+拉格朗日插值+主席树)


    题面

    (nleq 10^{12},kleq 100)

    题解

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

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

    大概是这样,设(g(n,j))表示(1)(n)之间的所有满足(i)是质数或者(i)的最小质因子大于(p_j)的所有(f(i))之和,我们根据递归地来求解(g),设一个阈值(L=sqrt{n}),当(nleq L)的时候,用主席树优化,能做到每一次查询只有(O(logsqrt{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=leftlfloorfrac{n}{m} ight floor),容易发现这里的(xleq sqrt{n}),想想整除分块,发现这里的(x)就是最大的满足(leftlfloorfrac{n}{x} ight floor=m)的数,那么每一个(x)都和(m)是一一对应的,我们就可以把它给存起来

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

    //minamoto
    #include<cstdio>
    #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 fd(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
    ",G(n,tot-1));
    	return 0;
    }
    
  • 相关阅读:
    软工实践个人总结
    第03组 每周小结 (3/3)
    第03组 每周小结 (2/3)
    第03组 每周小结(1/3)
    第03组 Beta冲刺 总结
    第03组 Beta冲刺 (5/5)
    第03组 Beta冲刺 (4/5)
    第03组 Beta冲刺 (3/5)
    第03组 Beta冲刺 (2/5)
    第03组 Beta冲刺 (1/5)
  • 原文地址:https://www.cnblogs.com/bztMinamoto/p/10415037.html
Copyright © 2020-2023  润新知