【JZOJ5224】【GDOI2018模拟7.12】C

Description

Data Constraint

Solution

首先必须讲讲自然数幂求和。

我们设Sk(n)=∑i=1nik

我们用第一类斯特林数来计算。

第一类Stirling数 s(p,k)

s(p,k)的一个的组合学解释是：将p个物体排成k个非空循环排列的方法数。

s(p,k)的递推公式： s(p,k)=(p-1)*s(p-1,k)+s(p-1,k-1) ,1<=k<=p-1

边界条件：s(p,0)=0 ,p>=1 s(p,p)=1 ,p>=0

递推关系的说明：

考虑第p个物品，p可以单独构成一个非空循环排列，这样前p-1种物品构成k-1个非空循环排列，方法数为s(p-1,k-1)。也可以前p-1种物品构成k个非空循环排列，而第p个物品插入第i个物品的左边，这有(p-1)*s(p-1,k)种方法。

由第一类斯特林数的定义可知Cpn=Ppnp!=S(p,p)np−S(p,p−1)np−1+……S(P,0)n0p!Cpn=Ppnp!=∑pi=0(−1)p−iS(p,i)nip!Sk(n)=∑j=1np!Cpj−∑i=0p−1(−1)p−iS(p,i)jiSp(n)=p!∑j=1nCpj−∑i=0p−1(−1)p−iS(p,i)∑j=1njiSp(n)=p!Cp+1n+1−∑i=0p−1(−1)p−iS(p,i)Si(n)Sp(n)=p!Cp+1n+1−∑i=0p−1(−1)p−iS(p,i)∑j=1njiSp(n)=(n−p+1)∗……∗(n+1)p+1−∑i=0p−1(−1)p−iS(p,i)Si(n)

由于我们知道S1(n)=n*(n+1)/2，所以我们可以在O(P *P)的时间内算出Sp(n)。

接下来回归正题。gcd(i,j)k=∑d=0ndk−∑i=1n/d∑j=1n/d[(i,j)==1]=∑d=0ndk−2∗∑i=1n/dϕ(i)−1前面∑nd=0dk用上面的方法求，后面的∑n/di=1ϕ(i)用杜教筛来求。不会杜教筛的请看欧拉函数之和。

时间复杂度O(N3/4)。由于本题模的地方较多，要尽量减少取模次数，卡卡常，否则会很难过……

Code

#include<iostream>
#include<cmath>
#include<cstring>
#include<cstdio>
#include<algorithm>
#define ll long long
using namespace std;
const ll maxn=8e6+5,mo=1e9+7,maxn1=maxn-5;
ll phi[maxn],h[maxn][2],sk[10],s[10][10];
int d[maxn],bz[maxn];
ll n,m,i,t,j,k,l,x,y,z,ans,p;
ll hash(ll x){
int t=x%maxn;
while (h[t][0] && h[t][0]!=x) t=(t+1)%maxn;
return t;
}
ll dg(ll n){
if (n<=maxn1) return phi
;
ll x=hash(n),t=n%mo,k,i=2,ans=0,y=(mo+1)/2;
if (h[x][0]) return h[x][1];h[x][0]=n;
ans=t*(t+1)%mo*y%mo;
while (i<=n){
t=n/i;k=n/(n/i);
ans=ans-(k-i+1)*dg(t)%mo;i=k+1;
}
h[x][1]=(ans%mo+mo)%mo;return h[x][1];
}
ll dg1(ll n){
ll x=n%mo,y=(mo+1)/2,p,i,j,k,t;
sk[1]=x*(x+1)%mo*y%mo;
for (p=2;p<=m;p++){
t=0;x=1;
for (j=n-p+1;j<=n+1;j++)
if (j%(p+1)==0 && !t) x=(j/(p+1))%mo*x%mo,t++;
else x=j%mo*x%mo;
for (j=0;j<p;j++){
k=((p-j)%2)?-1:1;
x=x-k*s[p][j]*sk[j]%mo;
}
sk[p]=(x%mo+mo)%mo;
}
return sk[m];
}
int main(){
//  freopen("data.in","r",stdin);
scanf("%lld%lld",&n,&m);phi[1]=1;
for (i=2;i<=maxn1;i++){
if (!bz[i]) d[++d[0]]=i,phi[i]=i-1;
for (j=1;j<=d[0];j++){
if (i*d[j]>maxn1) break;
bz[i*d[j]]=1;
if (i%d[j]==0){
phi[i*d[j]]=phi[i]*d[j];break;
}else phi[i*d[j]]=phi[i]*phi[d[j]];
}
}
for (i=1;i<=maxn1;i++)
phi[i]=(phi[i]+phi[i-1])%mo;
s[0][0]=1;
for (i=1;i<=m;i++){
for (j=1;j<i;j++)
s[i][j]=(s[i-1][j-1]+(i-1)*s[i-1][j])%mo;
s[i][i]=1;
}
i=1;
while (i<=n){
x=n/i;t=n/(n/i);y=2*dg(x)%mo-1;
ans=ans+(dg1(t)-dg1(i-1)+mo)%mo*y%mo;
i=t+1;
}
ans=ans%mo;
printf("%lld\n",ans);
}