EOJ Monthly 2019.11 E 数学题

题目链接

这题和多校以及今年网络赛的某题非常相似，算是比较套路的了。
先推f：
$\begin{aligned} f(n,k)&=\Sigma_{x_{1}=1}^{n}\Sigma_{x_{2}=1}^{n}\dots\Sigma_{x_{k}=1}^{n}[gcd(x_{1},x_{2},\dots,x_{k},n)=1] \\&=\Sigma_{x_{1}=1}^{n}\Sigma_{x_{2}=1}^{n}\dots\Sigma_{x_{k}=1}^{n}e(gcd(x_{1},x_{2},\dots,x_{k},n))\\&=\Sigma_{x_{1}=1}^{n}\Sigma_{x_{2}=1}^{n}\dots\Sigma_{x_{k}=1 }^{n}\Sigma_{d|gcd(x_{1},x_{2},\dots,x_{k},n)}\mu(d)\\&=\Sigma_{d|n}\mu(d)\Sigma_{x_{1}=1}^{\frac{n}{d}}\Sigma_{x_{2}=1}^{\frac{n}{d}}\dots\Sigma_{x_{k}=1 }^{\frac{n}{d}}1\\&=\Sigma_{d|n}\mu(d)(\frac{n}{d})^{k} \end{aligned}$f(n,k)​=Σx1​=1n​Σx2​=1n​…Σxk​=1n​[gcd(x1​,x2​,…,xk​,n)=1]=Σx1​=1n​Σx2​=1n​…Σxk​=1n​e(gcd(x1​,x2​,…,xk​,n))=Σx1​=1n​Σx2​=1n​…Σxk​=1n​Σd∣gcd(x1​,x2​,…,xk​,n)​μ(d)=Σd∣n​μ(d)Σx1​=1dn​​Σx2​=1dn​​…Σxk​=1dn​​1=Σd∣n​μ(d)(dn​)k​
所以：
$\begin{aligned} \Sigma_{i=1}^{n}f(i,k)&=\Sigma_{i=1}^{n}\Sigma_{d|n}\mu(d)(\frac{n}{d})^{k}\\&=\Sigma_{d=1}^{n}\mu(d)\Sigma_{i=1}^{\lfloor n/d \rfloor}i^{k} \end{aligned}$Σi=1n​f(i,k)​=Σi=1n​Σd∣n​μ(d)(dn​)k=Σd=1n​μ(d)Σi=1⌊n/d⌋​ik​
得到上式就好做了：后边的和式使用拉格朗日插值法计算前 $n$n 个数的 $k$k 次幂和，前边的和式使用 数论分块 以及 杜教筛 就可以完成这个式子的计算，复杂度大约为 $O(n^{\frac{1}{2}}k)$O(n21​k) 。

有关杜教筛莫比乌斯函数前缀和：
令
$S(n) = \Sigma_{i=1}^{n}\mu(i)$S(n)=Σi=1n​μ(i)则有
$\begin{aligned} 1&=\Sigma_{i=1}^{n}e(i)\\&= \Sigma_{i=1}^{n}\Sigma_{d|i}\mu(\frac{i}{d})I(d)\\&=\Sigma_{d=1}^{n}I(d)\Sigma_{i=1}^{\lfloor n/d \rfloor}\mu(i)\\&=\Sigma_{d=1}^{n}S(\lfloor n/d \rfloor)\\&=S(n)+\Sigma_{d=2}^{n}S(\lfloor n/d \rfloor) \end{aligned}$1​=Σi=1n​e(i)=Σi=1n​Σd∣i​μ(di​)I(d)=Σd=1n​I(d)Σi=1⌊n/d⌋​μ(i)=Σd=1n​S(⌊n/d⌋)=S(n)+Σd=2n​S(⌊n/d⌋)​即
$S(n)=1-\Sigma_{d=2}^{n}S(\lfloor n/d \rfloor)$S(n)=1−Σd=2n​S(⌊n/d⌋)

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<unordered_map>
#include<map>
using namespace std;
typedef long long ll;
const ll mod=998244353;
const int N=1005;
const int M=2000005;
ll n,k,mu[M]={},res[N]={};
ll coef[N]={},pref[N]={},suff[N]={};
unordered_map<ll,ll> mp;
//map<ll,ll> mp;
inline ll pow_mod(ll a,ll b) {
    ll ret=1;
    while(b) {
        if(b&1) ret=ret*a%mod;
        b>>=1;
        a=a*a%mod;
    }
    return ret;
}
inline ll mod_inv(ll a) {
    return pow_mod(a,mod-2);
}
inline void init(ll n,ll k) {
    res[0]=0;
    for(int i=1;i<=k+1;++i) {
        res[i]=res[i-1]+pow_mod(i,k);
        res[i]%=mod;
    }
    // calc coef[0]
    coef[0]=1;
    for(int i=-1;i>=-k-1;--i) {
        coef[0]*=(i+mod);
        coef[0]%=mod;
    }
    coef[0]=mod_inv(coef[0]);
    for(int i=1;i<=k+1;++i) {
        coef[i]=coef[i-1]*mod_inv(i)%mod;
        coef[i]*=(-k-2+i+mod);
        coef[i]%=mod;
    }
}
inline ll calc(ll n,ll k) {
    ll ret=0,tmp=1;
    if(n<=k+1) return res[n];
    if(n>=mod) {
        for(int i=0;i<=k+1;++i) {
            tmp=1;
            for(int j=0;j<=k+1;++j) {
                if(i!=j) {
                    tmp*=(n-j+mod)%mod;
                    tmp%=mod;
                }
            }
            ret+=tmp*coef[i]%mod*res[i]%mod;
            ret%=mod;
        }
        return ret;
    }
    else {
        pref[0]=n;
        //cout<<n<<" "<<k<<endl;
        for(int i=1;i<=k+1;++i) {
            pref[i]=pref[i-1]*((n-i+mod)%mod)%mod;
          //  cout<<pref[i]<<endl;
        }
        suff[k+1]=(n-k-1+mod)%mod;
        for(int i=k;i>=0;--i) {
            suff[i]=suff[i+1]*((n-i+mod)%mod)%mod;
            //cout<<suff[i]<<endl;
        }
        for(int i=0;i<=k+1;++i) {
            //cout<<coef[i]<<" "<<res[i]<<endl;
            if(!i) {
                ret+=suff[1]*coef[0]%mod*res[0]%mod;
            }
            else if(i==k+1) {
                ret+=pref[k]*coef[k+1]%mod*res[k+1]%mod;
            }
            else {
                ret+=suff[i+1]*pref[i-1]%mod*coef[i]%mod*res[i]%mod;
            }
            ret%=mod;
            //cout<<ret<<endl;
        }
    }
    return ret;
}
bool isprime[M]={};
ll prime[M]={},tot=0;
inline void getMu() {
    isprime[0]=isprime[1]=1;
    mu[0]=0; mu[1]=1;
    for(int i=2;i<M;++i) {
        if(!isprime[i]) {
            prime[tot++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<tot&&prime[j]*ll(i)<ll(M);++j) {
            isprime[prime[j]*i]=1;
            if(i%prime[j]==0) {
                mu[i*prime[j]]=0;
                break;
            }
            else {
                mu[i*prime[j]]=mu[i]*mu[prime[j]];
            }
        }
    }
    for(int i=1;i<M;++i) {
        mu[i]=(mu[i-1]+mu[i]+mod)%mod;
    }
}
inline ll getPre(ll n) {
    if(n<M) return mu[n];
    if(mp.find(n)!=mp.end()) {
        return mp[n];
    }
    ll r,ret=0;
    for(ll d=2;d<=n;++d) {
        r=n/(n/d);
        ret+=getPre(n/d)*(r-d+1)%mod;
        ret%=mod;
        d=r;
    }
    return mp[n]=(1-ret+mod)%mod;
}
inline ll check(ll n,ll k) {
    ll ret=0;
    for(int i=1;i<=n;++i) {
        ret+=pow_mod(i,k);
        ret%=mod;
    }
    return ret;
}
int main(void)
{
    getMu();
    scanf("%lld%lld",&n,&k);
    init(n,k);
    ll r=0,ans=0;
    for(ll d=1;d<=n;++d) {
        r=n/(n/d);
        ans+=(getPre(r)-getPre(d-1)+mod)%mod*calc(n/d,k)%mod;
        ans%=mod;
        d=r;
    }
    printf("%lld\n",ans);
    return 0;
}