【JZOJ 5683】【GDSOI2018模拟4.22】Prime

Description

求
$$\sum x^k$$
其中，x为n以内的质数，答案对1e9+7取模
$n\le 10^{12}$

Solution

用Min_25筛，拉格朗日差值算自然数幂和，
当然这是不够的，还要加上数据结构的优化，

对于递归的版本，每次我们要询问的是S(n,j)，在2~n中，每个x的$x^k$的和，满足x为质数或x最小质因子大于$p_j$，考虑当n小于阈值L时预处理，
对于L以内数，我们可以通过线筛筛出它的最小质因子，设x的最小质因子为$s{m}_x$，
那么上面的问题就是：询问2~n中，${\sum }_{i=2}^x{i}^k\ast [s{m}_i>p_j]$，
这个用主席树做即可，

复杂度：
可以发现这个过程和杜教筛很像，只是查询预处理的复杂度变成了log，
上界显然是$O(n^{\frac{2}{3}}+\sqrt{n}log(n))$，但远远达不到

如果有人知道更紧的上界麻烦留言，并不会QwQ

Code

#include <cstdio>
#include <algorithm>
#define fo(i,a,b) for(int i=a;i<=b;++i)
#define fod(i,a,b) for(int i=a;i>=b;--i)
#define min(q,w) ((q)>(w)?(w):(q))
#define max(q,w) ((q)<(w)?(w):(q))
using namespace std;
typedef long long LL;
const int N=2100500,mo=1e9+7;
int read(int &n)
{
	char ch=' ';int q=0,w=1;
	for(;(ch!='-')&&((ch<'0')||(ch>'9'));ch=getchar());
	if(ch=='-')w=-1,ch=getchar();
	for(;ch>='0' && ch<='9';ch=getchar())q=q*10+ch-48;n=q*w;return n;
}
int m;
LL n,ans;
int sum[N+1];
int pr[N/3+5],pr1[N+5],pri[N+5];
LL prp[N/3+5];
bool prz[N+5];
int jc[200],jcn[200],Fs[N+5];
int d[200];
struct SGT
{
	int l,r,v;
}b[N*18];
int root[N+1],b0;
void change(int l,int r,int &e,int l1,int l2)
{
	b[++b0]=b[e];((b[e=b0].v+=l2)>=mo?b[e].v-=mo:0);
	if(l==r)return;
	int t=(l+r)>>1;
	((l1<=t)?change(l,t,b[e].l,l1,l2):change(t+1,r,b[e].r,l1,l2));
}
LL find(int l,int r,int e,int l1)
{
	if(!e)return 0;
	if(l1<=l)return b[e].v;
	int t=(l+r)>>1;
	LL ans=find(t+1,r,b[e].r,l1);
	return((l1<=t?(ans+=find(l,t,b[e].l,l1)):ans)>=mo?ans-mo:ans);
}
int Fs1[N];
LL Gsum(LL n)
{
	if(n<=N)return Fs[n];
	int n1=(::n)/n;
	if(Fs1[n1])return Fs1[n1];
	n=n%mo;
	d[m+3]=1;fod(i,m+2,2)d[i]=d[i+1]*(n-(LL)i)%mo;
	LL ans=0,t=1;
	fo(i,1,m+2)
	{
		ans=(ans+t*d[i+1]%mo*(LL)jcn[i-1]%mo*(LL)jcn[m-i+2]%mo*(LL)Fs[i]%mo*(((m-i+2)&1)?-1LL:1LL))%mo;
		t=t*(n-(LL)i)%mo;
	}
	return Fs1[n1]=ans;
}
LL ksm(LL q,int w)
{
	LL ans=1;
	for(;w;w>>=1,q=q*q%mo)if(w&1)ans=ans*q%mo;
	return ans;
}
void Pre(int n)
{
	fo(i,2,n)
	{
		if(!prz[i])pr[++pr[0]]=i,Fs[i]=sum[i]=pr1[pr[0]]=ksm(i,m),prp[pr[0]]=(LL)i*i;
		fo(j,1,pr[0])
		{
			int t=i*pr[j];
			if(t>n)break;
			prz[t]=1;
			Fs[t]=(LL)Fs[i]*Fs[pr[j]]%mo;
			pri[t]=j;
			if(i%pr[j]==0)break;
		}
	}
	fo(i,2,n)
	{
		root[i]=root[i-1];
		if(pri[i])change(1,pr[0],root[i],pri[i],Fs[i]);
	}
	jc[0]=jc[1]=1;fo(i,2,m+10)jc[i]=jc[i-1]*(LL)i%mo;
	jcn[m+10]=ksm(jc[m+10],mo-2);fod(i,m+9,0)jcn[i]=jcn[i+1]*(i+1LL)%mo;
	fo(i,2,n)
	{
		((sum[i]+=sum[i-1])>=mo?sum[i]-=mo:0);
		((Fs[i]+=Fs[i-1])>=mo?Fs[i]-=mo:0);
	}
}
LL Gg(LL n,int m)
{
	if(n<N&&n<(LL)prp[m+1])return sum[n];
	if(n<=N)return (find(1,pr[0],root[n],m+1)+sum[n])%mo;
	if(!m)return Gsum(n);
	for(;n<(LL)prp[m];--m);
	LL ans=0;
	for(;m;--m)ans=(ans-(Gg(n/pr[m],m-1)-sum[pr[m]-1])*pr1[m])%mo;
	return (ans+Gsum(n))%mo;
}
int main()
{
	int q,w;
	scanf("%lld%d",&n,&m);
	Pre(N);
	printf("%lld\n",(Gg(n,pr[0]-1)+mo)%mo);
	return 0;
}