【jzoj5224】【GDOI2018模拟7.12】【C】【自然数幂和】【杜教筛】

题目大意

解题思路

$ans=\sum_{i=1}^n\sum_{j=1}^n(i,j)^k，k是给定的常数$

$ans=\sum_{d=1}^nd^k\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}[(i,j)==1]$

$=\sum_{d=1}^nd^k(2\sum_{j=1}^{n/d}\phi(j)-1)$

设$s(n)=\sum_{i=1}^n\phi(i)$

$ans==\sum_{d=1}^nd^k(2*s(n/d)-1)$

考虑分快处理

那么$s(n)=n(n+1)/2-\sum_{i=2}^ns(n/i)$

预处理出前一部分，记忆化后一部分即可。

对于前面的一部分使用自然数幂和求解。

$\sum_{i=0}^ni^t=\sum_{i=0}^tS(t,i)(n+1)^{（i+1(此处应有下划线)）}/(i+1)$

下划线处表示降幂

这样我们就可以快速地求出答案，虽然复杂度是对的然而我还是超时了两个点，表示常数卡得不行。

code

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define LF double
#define LL long long
#define ULL unsigned LL
#define fo(i,j,k) for(LL i=j;i<=k;i++)
#define fd(i,j,k) for(LL i=j;i>=k;i--)
#define fr(i,j) for(LL i=begin[j];i;i=next[i])
using namespace std;
LL max(LL x,LL y){return (x>y)?x:y;}
LL min(LL x,LL y){return (x<y)?x:y;}
LL const mn=1e6+9,si=1e6,mo=1e9+7;
LL lim=1e6,n,K,tag[mn],ss[mn],ans,s[mn],hs[mn],f[mn],g[10][10];
LL pow(LL x,LL y){
    LL z=1;
    while(y){
        if(y&1)z=z*x%mo;
        x=x*x%mo;
        y>>=1;
    }
    return z;
}
LL hash(LL x){
    int y=x%si;
    while(hs[y]&&(hs[y]!=x))y=(y+1==si)?0:y+1;
    return y;
}
LL mode(LL x){
    while(x>=mo)x-=mo;
    return x;
}
LL S(LL n){
    if(n<=lim)return s[n];
    LL tmp=hash(n);
    if(hs[tmp])return f[tmp];
    hs[tmp]=n;
    LL tm2=((n&1)==0)?mode(n/2)*mode(n+1)%mo:
        mode((n+1)/2)*mode(n%mo)%mo;
    for(LL i=2,j;i<=n;i=j+1){
        LL tmp=n/i;j=n/tmp;
        tm2=(tm2-mode(j-i+1)*S(tmp))%mo;
    }
    return f[tmp]=tm2;
}
LL calc(LL n,LL t){
    LL ans=0;
    fo(i,0,t){
        LL tmp=g[t][i];
        fd(j,n+1,n+1-i)tmp=tmp*(j%mo)%mo;
        ans=(ans+tmp*pow(i+1,mo-2))%mo;
    }
    return ans;
}
int main(){
    freopen("d.in","r",stdin);
    freopen("d.out","w",stdout);
    scanf("%lld%d",&n,&K);s[1]=1;
    g[0][0]=1;
    fo(i,1,5)fo(j,1,5)g[i][j]=g[i-1][j-1]+j*g[i-1][j];
    lim=1e6;
    fo(i,2,lim){
        if(!tag[i])ss[++ss[0]]=i,s[i]=i-1;
        fo(j,1,ss[0]){
            if(i*ss[j]>lim)break;
            tag[i*ss[j]]=1;
            if(i%ss[j]==0){s[i*ss[j]]=s[i]*ss[j];break;}
            else s[i*ss[j]]=s[i]*(ss[j]-1);
        }
        s[i]=(s[i]+s[i-1])%mo;
    }
    LL ans=0;
    for(LL i=1,j;i<=n;i=j+1){
        LL tmp=n/i;j=n/tmp;
        ans=(ans+(2*S(tmp)-1)*(calc(j,K)-calc(i-1,K)))%mo;
    }
    printf("%lld",(ans+mo)%mo);
    return 0;
}