【JZOJ5224】【GDOI2018模拟7.12】C

Description

Data Constraint

Solution

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

$$Sk(n)=\sum_{i=1}^ni^k$$
我们用第一类斯特林数来计算。
第一类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)种方法。
由第一类斯特林数的定义可知 $$C^p_n={P^p_n\over p!}={S(p,p)n^p-S(p,p-1)n^{p-1}+……S(P,0)n^0\over p!}$$ $$C^p_n={P^p_n\over p!}={\sum_{i=0}^p(-1)^{p-i}S(p,i)n^i\over p!}$$ $$Sk(n)=\sum_{j=1}^np!C_j^p-\sum_{i=0}^{p-1}(-1)^{p-i}S(p,i)j^i$$ $$Sp(n)=p!\sum_{j=1}^nC_j^p-\sum_{i=0}^{p-1}(-1)^{p-i}S(p,i)\sum_{j=1}^nj^i$$ $$Sp(n)=p!C_{n+1}^{p+1}-\sum_{i=0}^{p-1}(-1)^{p-i}S(p,i)Si(n)$$ $$Sp(n)=p!C_{n+1}^{p+1}-\sum_{i=0}^{p-1}(-1)^{p-i}S(p,i)\sum_{j=1}^nj^i$$ $$Sp(n)={(n-p+1)*……*(n+1)\over p+1}-\sum_{i=0}^{p-1}(-1)^{p-i}S(p,i)Si(n)$$
由于我们知道S1(n)=n*(n+1)/2，所以我们可以在O(P *P)的时间内算出Sp(n)。
接下来回归正题。 $$gcd(i,j)^k=\sum_{d=0}^{n}d^k-\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}[(i,j)==1]$$ $$=\sum_{d=0}^{n}d^k-2*\sum_{i=1}^{n/d}\phi(i)-1$$ 前面 $\sum_{d=0}^{n}d^k$用上面的方法求，后面的 $\sum_{i=1}^{n/d}\phi(i)$用杜教筛来求。不会杜教筛的请看 欧拉函数之和。
时间复杂度O( $N^{3/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[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);
}