EOJ Monthly 2019.11 E 数学题

EOJ Monthly 2019.11 E 数学题

题目链接

这题和多校以及今年网络赛的某题非常相似，算是比较套路的了。
先推f：
$$\begin{array}{rl}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{array}$$
所以：
$$\begin{array}{rl}{\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{array}$$
得到上式就好做了：后边的和式使用拉格朗日插值法计算前 $n$ 个数的 $k$ 次幂和，前边的和式使用 数论分块 以及 杜教筛 就可以完成这个式子的计算，复杂度大约为 $O(n^{\frac{1}{2}}k)$ 。

有关杜教筛莫比乌斯函数前缀和：
令
$$S(n)={\Sigma }_{i=1}^n\mu (i)$$则有
$$\begin{array}{rl}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{array}$$即
$$S(n)=1-{\Sigma }_{d=2}^{n}S(\lfloor n/d\rfloor )$$

代码

#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;
}