EOJ Monthly 2019.11 E. 数学题（莫比乌斯反演）

题意

输入\(n,k\),求\(\sum_{n=1}^{N}\sum_{a_1}^n...\sum_{a_k}^n[gcd(a_1,...a_k,n)==1]\)

思路

设\(f(d)=\sum_{a_1}^N...\sum_{a_k}^N[gcd(a_1,...a_k,N)==d] \\ F(n)=\sum_{n|d}f(d)=\lfloor\frac{N}{n}\rfloor^k \ (n|N)\\ f(n)=\sum_{n|d}\mu(\frac{d}{n})F(d) \\ f(1)=\sum_{d=1}^N\mu(d)F(d)=\sum_{d|N}\mu(d)(\frac{N}{d})^k \\ ans=\sum_{n=1}^{N}\sum_{d|n}\mu(d)(\frac{n}{d})^k=(枚举d)\sum_{d=1}^{N}\mu(d)\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}i^k\)
用分块+杜教筛求该式子，幂方可用拉格朗日插值法求也可用伯努利数去求

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int maxx = 1e9+10;
const int N = 1e6+10;
const int K = 1010;
const int mod = 998244353;
LL inv[K],b[K],c[K][K],tmp[K];
LL phi[N],vis[N],prime[N],mu[N],sum[N];
map<LL,LL>w;
int n,k;
void getmu()
{
    int cnt=0;
    phi[1]=mu[1]=1;
    for(int i=2;i<N;i++)
    {
        if(!vis[i])prime[++cnt]=i,mu[i]=-1,phi[i]=i-1;
        for(int j=1;j<=cnt&&prime[j]*i<N;j++)
        {
            vis[i*prime[j]]=1;
            if(i%prime[j]==0)
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            else mu[i*prime[j]]=-mu[i],phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
    }
    for(int i=1;i<N;i++)sum[i]=sum[i-1]+mu[i];
}
int dfsmu(int x)
{
    if(x<=1e6)return sum[x];
    if(w[x])return w[x];
    int ans=1;
    for(int l=2,r;l>=0&&l<=x;l=r+1)
    {
        r=x/(x/l);
        ans-=(r-l+1)*dfsmu(x/l);
    }
    return w[x]=ans;
}
void getb()
{
    for(int i=0;i<K;i++)
    {
        c[i][0]=1;
        for(int j=1;j<=i;j++)
            c[i][j]=(c[i-1][j]%mod+c[i-1][j-1]%mod)%mod;
    }
    inv[1]=1;
    for(int i=2;i<K;i++)
        inv[i]=(mod-mod/i)*inv[mod%i]%mod;
    b[0]=1;
    for(int i=1;i<K;i++)
    {
        LL ans=0;
        if(i==K-1)break;
        for(int j=0;j<i;j++)
            ans=(ans+c[i+1][j]*b[j])%mod;
        ans*=-inv[i+1];
        ans=(ans%mod+mod)%mod;
        b[i]=ans;
    }
}
LL work(int t)
{
    tmp[0]=1;
    for(int i=1;i<=k+1;i++)
        tmp[i]=tmp[i-1]*(t+1)%mod;
    LL ans=inv[k+1];
    LL sum=0;
    for(int i=1;i<=k+1;i++)
    {
        sum+=c[k+1][i]*tmp[i]%mod*b[k+1-i]%mod;
        sum%=mod;
    }
    ans*=sum;
    ans%=mod;
    return ans;
}
int main()
{
    getb();
    getmu();
    scanf("%d%d",&n,&k);
    LL ans=0;
    for(int l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        ans=(ans+work(n/l)*(dfsmu(r)-dfsmu(l-1))%mod)%mod;
    }
    ans=(ans%mod+mod)%mod;
    printf("%lld\n",ans);
}