C. 在表格里造序列

题意

对于每一对满足 \(1\le i,j\le n\) 的 \((i,j)\)，计算有多少个长度为 \(m\) 的序列，权值在 \([1,n]\) 之间且 \(\gcd(a_1,a_2,...,a_m)=\gcd(i,j)\)。答案对 \(998244353\) 取模。

---

思路

方法：莫比乌斯反演+杜教筛

不会莫比乌斯反演？ 出门右转：OI-wiki。
不会杜教筛？ 出门右转：OI-wiki。

现在进入正文部分。

首先，我们考虑他的子问题：计算有多少个长度为 \(m\) 的序列，权值在 \([1,n]\) 之间且 \(\gcd(a_1,a_2,...,a_m)=k\)。

显然的，\(a_i\) 一定是 \(k\) 的倍数。则我们给所有数同时除 \(k\)，现在就权值就变成了 \(\left[1,\left\lfloor\frac nk\right\rfloor\right]\)，限制也变成了 \(\gcd(a_1,a_2,...,a_m)=1\)。

我们设 \(f_i\) 表示 \(\gcd(a_1,a_2,...,a_m)=i\) 的方案数，\(g_i\) 表示仅仅满足 \(i|a_j\) 的方案数。这也意味着 \(g_i\) 计算的是 \(i|\gcd(a_1,a_2,...,a_m)\) 的方案数。

那么显然的，\(g_k=\left\lfloor\frac nk\right\rfloor^m\)。而且，\(f\) 和 \(g\) 还满足以下关系式：

\[g_k=\sum_{k|d}f_d \]

那么根据莫比乌斯反演的式子，则有：

\[f_k=\sum_{k|d}\mu\left(\frac dk\right)g(d) \]

我们令 \(d=ik\)，则：

\[f_k=\sum_{i=1}^{\left\lfloor n/k\right\rfloor}\mu(i)g(ik)=\sum_{i=1}^{\left\lfloor n/k\right\rfloor}\mu(i)\left\lfloor\frac n{ik}\right\rfloor^m \]

现在我们知道了子问题的答案，那么总问题的呢？

答案就很显然了！就是这个式子：

\[\sum_{i=1}^n\sum_{j=1}^nf_{\gcd(i,j)} \]

当然，这个式子肯定是过不了的。于是我们考虑推式子。

\[\begin{align*} &\sum_{i=1}^n\sum_{j=1}^nf_{\gcd(i,j)}\\ =&\sum_{k=1}^n\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=k]f_k\\ =&\sum_{k=1}^nf_k\sum_{i=1}^{\left\lfloor n/k\right\rfloor}\sum_{j=1}^{\left\lfloor n/k\right\rfloor}[\gcd(i,j)=1]\\ \end{align*} \]

现在又要用到狄利克雷卷积的基本式：\(\epsilon=\mu*1\)。也就是说：

\[\sum_{d|n}\mu(d)=[n=1] \]

所以我们的式子又变成了这个样子：

\[\begin{align*} &\sum_{k=1}^nf_k\sum_{i=1}^{\left\lfloor n/k\right\rfloor}\sum_{j=1}^{\left\lfloor n/k\right\rfloor}[\gcd(i,j)=1]\\ =&\sum_{k=1}^nf_k\sum_{i=1}^{\left\lfloor n/k\right\rfloor}\sum_{j=1}^{\left\lfloor n/k\right\rfloor}\sum_{d|\gcd(i,j)}\mu(d)\\ =&\sum_{k=1}^nf_k\sum_{i=1}^{\left\lfloor n/k\right\rfloor}\sum_{j=1}^{\left\lfloor n/k\right\rfloor}\sum_{d|i,d|j}\mu(d)\\ =&\sum_{k=1}^nf_k\sum_{d=1}^{\left\lfloor n/k\right\rfloor}\sum_{i=1}^{\left\lfloor n/kd\right\rfloor}\sum_{j=1}^{\left\lfloor n/kd\right\rfloor}\mu(d)\\ =&\sum_{k=1}^nf_k\sum_{d=1}^{\left\lfloor n/k\right\rfloor}\mu(d)\left\lfloor\frac{n}{kd}\right\rfloor^2\\ =&\sum_{k=1}^n\left(\sum_{i=1}^{\left\lfloor n/k\right\rfloor}\mu(i)\left\lfloor\frac n{ik}\right\rfloor^m\right)\left(\sum_{d=1}^{\left\lfloor n/k\right\rfloor}\mu(d)\left\lfloor\frac{n}{kd}\right\rfloor^2\right)\\ \end{align*} \]

发现了吗？ 两个式子竟长的如此相似！

于是我们考虑里面的式子可以用整除分块做。而观察到每一个除法都带有 \(\left\lfloor\frac nk\right\rfloor\) 的因式，于是也可以对 \(k\) 做整除分块。

因为 \(n\) 有 \(10^9\)，于是要用杜教筛算 \(\mu\) 函数。

复杂度 \(O(n^{0.75})\)。

---

代码

#include <bits/stdc++.h>
#pragma GCC optimize(2)

using namespace std;

#define map unordered_map
#define int long long

#define fileio(x,y) freopen(x,"r",stdin),freopen(y,"w",stdout)
#define tup tuple<int,int,int>
#define pii pair<int,int>
#define bit(x) bitset<x>
#define pb emplace_back
#define mt make_tuple 
#define mp make_pair

const int   mod=998244353;
const int   maxn=2e6+10;

map<int,int>    mu,pp;
bit(maxn)       vis;
int             summu[maxn];
int             pri[maxn];
int             n,m,cnt;

int power(int a,int b,int p) { int tar=1; for(; b; b>>=1,a=1ll*a*a%p) if(b&1) tar=1ll*tar*a%p; return tar; }

int calcMu(int x)
{
    if(x<maxn) return summu[x];
    if(mu[x]) return mu[x];
    int val=1;
    for(int l=2,r; l<=x; l=r+1)
    {
        r=(x/(x/l));
        (val-=calcMu(x/l)*(r-l+1)%mod)%=mod;
    }
    return (mu[x]=val);
}

void init()
{
    summu[1]=1;
    for(int i=2; i<maxn; i++)
    {
        if(!vis[i]) { pri[++cnt]=i; summu[i]=-1; }
        for(int j=1; j<=cnt&&i*pri[j]<maxn; j++)
        {
            vis[i*pri[j]]=true;
            if(i%pri[j]==0) { summu[i*pri[j]]=0; break; }
            summu[i*pri[j]]=-summu[i];
        }
    }
    for(int i=1; i<maxn; i++) (summu[i]+=summu[i-1])%=mod;
    return;
}

void work()
{
    /* Code */
    cin>>n>>m;
    // n=1e9,m=1e18;
    int tar=0; m%=(mod-1);
    for(int l=1,r; l<=n; l=r+1) { r=(n/(n/l)); pp[n/l]=power(n/l,m,mod); }
    for(int l=1,r; l<=n; l=r+1)
    {
        r=(n/(n/l));
        // initialization
        int temp=n/l,tar1=0,tar2=0;
        for(int p=1,q; p<=temp; p=q+1)
        {
            q=(temp/(temp/p));
            int qmu=calcMu(q)-calcMu(p-1);
            (tar1+=pp[temp/p]*qmu%mod)%=mod;
            (tar2+=qmu*(temp/p)%mod*(temp/p)%mod)%=mod;
        }
        tar+=(r-l+1)*tar1%mod*tar2%mod;
    }
    cout<<(tar%mod+mod)%mod<<'\n';
    return;
}

signed main()
{
    // fileio(".in",".out");
    ios::sync_with_stdio(false);
    cin.tie(0);
    int t,stTime=clock();
    init();
    t=1;
    while(t--) work();
    // cerr<<"Time : "<<(double)(clock()-stTime)/CLOCKS_PER_SEC<<'\n';
    return 0;
} // Texas yyds !!!