【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)ji

Sp(n)=p!∑j=1nCpj−∑i=0p−1(−1)p−iS(p,i)∑j=1nji

Sp(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=1nji

Sp(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 longusing 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[n];    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);}