Description
求
∑
x
k
\sum x^k
∑xk
其中,x为n以内的质数,答案对1e9+7取模
n
≤
1
0
12
n\leq 10^{12}
n≤1012
Solution
用Min_25筛,拉格朗日差值算自然数幂和,
当然这是不够的,还要加上数据结构的优化,
对于递归的版本,每次我们要询问的是S(n,j),在2~n中,每个x的
x
k
x^k
xk的和,满足x为质数或x最小质因子大于
p
j
p_j
pj,考虑当n小于阈值L时预处理,
对于L以内数,我们可以通过线筛筛出它的最小质因子,设x的最小质因子为
s
m
x
sm_x
smx,
那么上面的问题就是:询问2~n中,
∑
i
=
2
x
i
k
∗
[
s
m
i
>
p
j
]
\sum_{i=2}^xi^k*[sm_i>p_j]
∑i=2xik∗[smi>pj],
这个用主席树做即可,
复杂度:
可以发现这个过程和杜教筛很像,只是查询预处理的复杂度变成了log,
上界显然是
O
(
n
2
3
+
n
l
o
g
(
n
)
)
O(n^{\frac{2}{3}}+\sqrt{n}log(n))
O(n32+nlog(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;
}