题解 P6384 【『MdOI R2』Quo Vadis】
Karry5307
·
·
题解
题意
题面写的非常清晰，不解释，摆数据范围太麻烦
题解
永远滴神 EA 的题解可能含有一点错误，写一个至少看上去比较对的官方题解。
日常给大家推荐 M2U 的歌呢/cy
首先直接证明结论，消元后的 A_{i,j} 可以表示成
A^\prime_{i,j}=
\begin{cases}
ij\varphi(i),&i\mid j
\\
0,&i\nmid j
\end{cases}
考虑利用一下高斯消元的性质（线性代数常识）提出每行每列的 i 和 j，构造新的矩阵 B_{i,j}=\gcd(i,j)，那么消元后的 B_{i,j} 为
B^\prime_{i,j}=
\begin{cases}
\varphi(i),&i\mid j
\\
0,&i\nmid j
\end{cases}
这个是 200 年前某篇论文的研究成果，这里简要的给大家证明一下。
考虑数学归纳法。
由于上述结论在 1 到 n-1 行成立，考虑第 n 行，注意到只有 n 的约数行才会给 n 带来贡献，于是我们得到
B^\prime_{n,k}=B_{n,k}-\sum\limits_{d\mid n,d<n}\varphi(d)
讨论一下，当 n\nmid k 的时候，B_{n,k}=\gcd(n,k)
此时注意到 \gcd(n,k)\mid n，所以 B_{\gcd(n,k),k} 一定会给 B_{n,k} 带来的贡献，并且系数为 1。
对于 \gcd(n,k) 前面的一些行数给 B^\prime_{n,k} 和 B^\prime_{\gcd(n,k),k} 带来相同的贡献，所以消元到第 \gcd(n,k)-1 行的时候有 B^\prime_{n,k}=B^\prime_{\gcd(n,k),k}。
于是第 \gcd(n,k) 行向第 n 行贡献完之后 B^\prime_{n,k} 就为 0。
否则考虑这样一个事实，每一行消元结束的时候都有 B^\prime_{n,k}=B^\prime_{n,n}
所以根据上式，得到
B^\prime_{n,k}=B^\prime_{n,n}=n-\sum\limits_{d\mid n,d<n}\varphi(d)=n+\varphi(n)-\sum\limits_{d\mid n}\varphi(d)
于是根据狄利克雷卷积的基本公式得到
B^\prime_{n,k}=\varphi(n)
同时该结论在第一行也成立，于是结论就对任何情况下都成立了。
然后我们来考虑每一个操作。
$1$ 操作前的 $3$ 操作是 P3768，也就是让我们求这个东西（为了复习莫比乌斯反演就把过程写在这里）
$$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}ij\gcd(i,j)$$
枚举 $\gcd$ 并且提出来得到
$$\sum\limits_{d=1}^{n}d^3\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}ij[\gcd(i,j)=1]$$
后面是套路反演
$$\sum\limits_{d=1}^{n}d^3\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}ij\sum\limits_{x\mid\gcd(i,j)}\mu(x)$$
提出来一下什么的
$$\sum\limits_{d=1}^{n}d^3\sum\limits_{x=1}^{\lfloor\frac{n}{d}\rfloor}x^2\mu(x)\left(\sum\limits_{i=1}^{\lfloor\frac{n}{xd}\rfloor}i\right)^2$$
我们设 $f(n)=\left(\sum\limits_{i=1}^{n}i\right)^2$，则有
$$\sum\limits_{d=1}^{n}d\sum\limits_{x=1}^{\lfloor\frac{n}{d}\rfloor}(xd)^2\mu(x)f(\lfloor\frac{n}{xd}\rfloor)$$
接下来记 $T=xd$，得到
$$\sum\limits_{T}^{n}T^2f(\lfloor\frac{n}{T}\rfloor)\sum\limits_{d\mid T}d\mu(\frac{T}{d})$$
右边是 $id$ 与 $\mu$ 的狄利克雷卷积，所以答案为
$$\sum\limits_{T}^{n}T^2f(\lfloor\frac{n}{T}\rfloor)\varphi(T)$$
这个可以整除分块和杜教筛搞定。
然后 $1$ 操作后的 $3$ 操作，考虑一列一列的搞，设 $g(n)$ 为第 $n$ 列的答案，得到
$$g(n)=n\sum\limits_{d\mid n}d\varphi(d)$$
设 $h(n)=\sum\limits_{d\mid n}d\varphi(d)$，对于质数 $p$ 和正整数 $k$ 我们有
$$h(p^k)=\sum\limits_{i=0}^{k}p^i\varphi(p^i)=\frac{p^{2k+1}+1}{p+1}$$
然后通过对 $n$ 质因数分解可以证明 $h(n)$ 是积性函数，于是线性筛即可。
行列式的话答案为 $(n!)^2\prod\limits_{i=1}^{n}\varphi(i)$。
### 代码
```cpp
#include<bits/stdc++.h>
using namespace std;
typedef int ll;
typedef long long int li;
const ll MAXN=5e6+51,MOD=998244353,INV6=166374059,INV4=748683265;
unordered_map<ll,ll>resP;
ll n,op,gauss,ptot;
li x,y;
ll prime[348521],phi[MAXN],np[MAXN],p[MAXN],low[MAXN],f[MAXN];
ll invPr[348521],prefixF[MAXN],prod[MAXN];
inline li read()
{
register li num=0,neg=1;
register char ch=getchar();
while(!isdigit(ch)&&ch!='-')
{
ch=getchar();
}
if(ch=='-')
{
neg=-1;
ch=getchar();
}
while(isdigit(ch))
{
num=(num<<3)+(num<<1)+(ch-'0');
ch=getchar();
}
return num*neg;
}
inline ll qpow(ll base,ll exponent)
{
ll res=1;
while(exponent)
{
if(exponent&1)
{
res=(li)res*base%MOD;
}
base=(li)base*base%MOD,exponent>>=1;
}
return res;
}
inline void sieve(ll limit)
{
ll pp;
phi[1]=f[1]=prod[0]=1;
for(register int i=2;i<=limit;i++)
{
if(!np[i])
{
low[i]=prime[++ptot]=i,phi[i]=i-1,f[i]=((li)i*i%MOD-i+1+MOD)%MOD;
invPr[ptot]=qpow(prime[ptot]+1,MOD-2);
}
for(register int j=1;j<=ptot&&prime[j]*i<=limit;j++)
{
np[prime[j]*i]=1;
if(!(i%prime[j]))
{
low[i*prime[j]]=low[i]*prime[j];
phi[i*prime[j]]=phi[i]*prime[j];
if(low[i]==i)
{
pp=i*prime[j];
f[pp]=(li)((li)pp*pp%MOD*prime[j]%MOD+1)*invPr[j]%MOD;
}
else
{
f[i*prime[j]]=(li)f[i/low[i]]*f[low[i]*prime[j]]%MOD;
}
break;
}
low[i*prime[j]]=prime[j],f[i*prime[j]]=(li)f[i]*f[prime[j]]%MOD;
phi[i*prime[j]]=phi[i]*phi[prime[j]];
}
}
for(register int i=1;i<=limit;i++)
{
p[i]=(p[i-1]+(li)i*i%MOD*phi[i]%MOD)%MOD;
prefixF[i]=(prefixF[i-1]+(li)i*f[i]%MOD)%MOD;
prod[i]=(li)prod[i-1]*i%MOD*i%MOD*phi[i]%MOD;
}
}
inline ll calc2(ll x)
{
return (li)x*(x+1)%MOD*(2*x+1)%MOD*INV6%MOD;
}
inline ll calc3(ll x)
{
return (li)x*(x+1)%MOD*x%MOD*(x+1)%MOD*INV4%MOD;
}
inline ll prefixP(ll num)
{
if(num<=5e6)
{
return p[num];
}
if(resP[num])
{
return resP[num];
}
ll res=calc3(num);
for(register int l=2,r;l<=num;l=r+1)
{
r=num/(num/l);
res=(res-(li)prefixP(num/l)*(calc2(r)-calc2(l-1)+MOD)%MOD+MOD)%MOD;
}
return resP[num]=res;
}
inline ll calc(ll num)
{
ll res=0;
for(register int l=1,r;l<=num;l=r+1)
{
r=num/(num/l);
res=(res+(li)calc3(num/l)*(prefixP(r)-prefixP(l-1)+MOD)%MOD)%MOD;
}
return res;
}
inline ll getPhi(ll x)
{
ll res=x;
for(register int i=2;i<=sqrt(x);i++)
{
if(x%i==0)
{
res=res/i*(i-1);
while(x%i==0)
{
x/=i;
}
}
}
if(x!=1)
{
res=res/x*(x-1);
}
return res;
}
int main()
{
sieve(5000010),n=read();
for(register int i=0;i<n;i++)
{
op=read();
if(op==1)
{
gauss=1;
}
if(op==2)
{
x=read(),y=read();
if(gauss)
{
printf("%lld\n",y%x?0:(li)x*y%MOD*getPhi(x)%MOD);
}
else
{
printf("%lld\n",(li)(x%MOD)*(y%MOD)%MOD*(__gcd(x,y)%MOD)%MOD);
}
}
if(op==3)
{
x=read();
if(gauss)
{
printf("%d\n",prefixF[x]);
}
else
{
printf("%d\n",calc(x));
}
}
if(op==4)
{
x=read();
printf("%d\n",prod[x]);
}
}
}
```