EOJ Monthly 2019.11 - E：数学题【反演+杜教筛+拉格朗日插值】

题目：

EOJ Monthly 2019.11 - E：数学题

题意：

咕咕咕

分析：

题目等价于求解：

$$\sum _{i=1}^n\sum _{a_1=1}^i\sum _{a_2=1}^i...\sum _{a_k=1}^i[gcd(a_1,a_2,...,a_n,i)==1]$$

$$\sum _{i=1}^n\sum _{d|i}\sum _{a_1=1}^{\lfloor \frac{i}{d}\rfloor }\sum _{a_2=1}^{\lfloor \frac{i}{d}\rfloor }...\sum _{a_k=1}^{\lfloor \frac{i}{d}\rfloor }\mu (d)$$$$\sum _{i=1}^n\sum _{d|i}\mu (d)(\lfloor \frac{i}{d}\rfloor {)}^k$$
令 $f(n)={\sum }_{d|n}\mu (d)(\lfloor \frac{n}{d}\rfloor {)}^k$，$g(n)=n^k$，观察到 $f(n)$ 为狄利克雷卷积展开形式，即：$f(n)=\mu \ast g(n)$，且为积性函数，题目所求转化为：$$S(n)=\sum _{i=1}^{n}f(i)$$
我们知道 $\mu \ast I=ϵ$，设 $h(n)=f\ast I(n)=ϵ\ast g(n)=g(n)$，有：
$$\sum _{i=1}^{n}h(i)=\sum _{i=1}^n\sum _{d|i}f(\frac{i}{d})I(d)=\sum _{d=1}^{n}I(d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }f(i)$$
提出等式右边的第一项：
$$\sum _{i=1}^{n}h(i)=I(1)S(n)+\sum _{d=2}^{n}I(d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }f(i)$$
恒等函数 $I(n)=1$，$S(n)$ 移到等式左边：
$$S(n)=\sum _{i=1}^{n}h(i)-\sum _{d=2}^n\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }f(i)$$$$S(n)=\sum _{i=1}^n{i}^k-\sum _{d=2}^{n}S(\lfloor \frac{n}{d}\rfloor )$$
$k$ 次幂前缀和用拉格朗日插值 $O(k)$ 计算得到，套用杜教筛 $O(n^{\frac{2}{3}}+k\ast \sqrt{n})$ 即可求解

代码：

#include <bits/stdc++.h>
 
#define x first
#define y second
#define pii pair<int,int>
#define sz(x) (int)(x).size()
#define Max(x,y) (x)>(y)?(x):(y)
#define Min(x,y) (x)<(y)?(x):(y)
#define all(x) (x).begin(),(x).end()
using namespace std;
typedef long long LL;
const int maxn = 2e6+56;
const int mod = 998244353;
LL qpow(LL a,LL x,LL mod){
    LL res = 1ll;
    while(x){
        if(x&1) res = res * a % mod;
        a = a * a % mod;
        x >>= 1;
    }
    return res;
}
LL per[maxn],suf[maxn],fac[maxn],a[maxn];
void init1(int n,int mod){
    LL res = 1;
    for(int i = 1;i <= n+2; ++i) res = res*i%mod;
    fac[n+2] = qpow(res,mod-2,mod);
    for(int i = n+1;i >= 0; --i) fac[i] = (i+1)*fac[i+1]%mod;
}
//给定最高次为n的多项式的n+1个点分别为：(0,f0),(1,f1),(2,f2)...(n,fn),求f(k)的值
LL Lagrange(LL f[],int n,int k){                   
    if(k <= n) return f[k];
    per[0] = suf[n+1] = 1;
    for(int i = 0;i <= n; ++i) per[i+1] = per[i]*(k-i)%mod;
    for(int i = n;i >= 0; --i) suf[i] = suf[i+1]*(k-i)%mod;
    LL fk = 0;
    for(int i = 0;i <= n; ++i){
        LL tep = f[i]*per[i]%mod*suf[i+1]%mod*fac[i]%mod*fac[n-i]%mod;
        if((n-i)&1) fk = (fk-tep+mod)%mod;
        else fk = (fk+tep)%mod; 
    }
    return fk;
}
int prime[maxn],vis[maxn],mul[maxn],f[maxn];
void init(int k){
    int tot = 0; mul[1] = 1; a[1] = 1;
    for(int i = 2;i < maxn; ++i){
        a[i] = qpow(i,k,mod);
        if(!vis[i]) prime[tot++] = i,mul[i] = -1;
        for(int j = 0;j<tot&&1ll*i*prime[j]<maxn; ++j){
            int v = i*prime[j]; vis[v] = 1;
            if(i%prime[j] == 0){
                mul[v] = 0; break;
            }
            mul[v] = -mul[i];
        }
    }
    for(int i = 1;i < maxn; ++i){
        for(int j = i;j < maxn;j += i){
            f[j] += (1ll*mul[i]*a[j/i]%mod+mod)%mod;
            if(f[j] >= mod) f[j] -= mod;
        }
        f[i] += f[i-1];
        if(f[i] >= mod) f[i] -= mod;
    } 
}
map<int,int> mp;
int solve(int n,int k){
    if(n < maxn) return f[n];
    if(mp.count(n)) return mp[n]; 
    int ans = 0;
    for(int i = 2,last;i <= n;i = last+1){
        int t = n/i; last = n/t;
        ans += 1ll*(last-i+1)*solve(t,k)%mod;
        if(ans >= mod) ans -= mod;
    }
    return mp[n] = (Lagrange(a,k+1,n)-ans+mod)%mod;
}
int main(){
    int n,k; scanf("%d %d",&n,&k);
    init(k); init1(k+5,mod);
    for(int i = 1;i <= k+1; ++i){
        a[i] += a[i-1];
        if(a[i] >= mod) a[i] -= mod;
    }
    printf("%d\n",solve(n,k));
    return 0;
}