UR#5 怎样跑得更快

给你整数 \(n,c,d\) ，现在有整数 \(x_1,…,x_n\) 和 \(b_1,…,b_n\)

满足：

 

\[\sum_{j=1}^n\gcd(i,j)^clcm(i,j)^dx_j≡b_i\pmod p \]

 

给出 \(b_1,…,b_n\)，请你解出 \(x_1,…,x_n\) 的值。

直接推式子，有

 

\[\begin{aligned} \sum_{j=1}^n\gcd(i,j)^clcm(i,j)^dx_j&=\sum_{j=1}^n\gcd(i,j)^{c-d}i^dj^dx_j\\ &=i^d\sum_{j=1}^n\gcd(i,j)^{c-d}j^dx_j\\ \end{aligned} \]

 

令 \(x_j=j^dx_j\)，这样最后再除掉一个 \(j^d\) 即可，所以现在得到了一个系数是和 \(\gcd\) 有关的方程组，令 \(F(x)=x^{c-d}\)，那么将方程组写成一个矩阵的形式，大概是下边这样的。

 

\[\begin{bmatrix} F(\gcd(1,1))& F(\gcd(1,2)) & F(\gcd(1,3)) & F(\gcd(1,4)) \\ F(\gcd(2,1)) & F(\gcd(2,2)) & F(\gcd(2,3)) & F(\gcd(2,4)) \\ F(\gcd(3,1)) & F(\gcd(3,2)) & F(\gcd(3,3)) & F(\gcd(3,4)) \\ F(\gcd(4,1)) & F(\gcd(4,2)) & F(\gcd(4,3)) & F(\gcd(4,4)) \end{bmatrix} \]

 

暴力做当然是直接高斯消元，不过注意到对于这种矩阵，一个事情是它可以直接消元消成一种比较好看的形式。

令 \(F(x)=\sum_{d|x}G(d)\)，那么原来的第 \(i\) 行第 \(j\) 列变成了 \(\sum_{d|i,d|j}G(d)\)，现在证明经过若干次初等行变换可以得到一个新矩阵，该矩阵的第 \(i\) 行第 \(j\) 列为 \([i|j]G(i)\) 。

使用归纳法证明，显然第一行不用任何变换就满足，对于其它行，只需要减去除了它以外的所有因子的行即可。

所以最后可以得到一个比较简单的矩阵。

 

\[\begin{bmatrix} G(\gcd(1,1))& G(\gcd(1,2)) & G(\gcd(1,3)) & G(\gcd(1,4)) \\ 0 & G(\gcd(2,2)) & 0 & G(\gcd(2,4)) \\ 0 & 0 & G(\gcd(3,3)) & 0 \\ 0 & 0 & 0 & G(\gcd(4,4)) \end{bmatrix} \]

 

对这个矩阵手动消一下元就能够得到每一个 \(x_j\) 的值。

判断无解和多组解的方法跟高斯消元一样。

#include<bits/stdc++.h>
using namespace std;
#define rint register int
#define ll long long
#define rll register long long
const int N=1e5+10;
const int p=998244353;
ll fpow(rll a,rll b){
	rll res=1;
	for(;b;b>>=1,a=a*a%p)
		if(b&1)res=res*a%p;
	return res;
}
ll b[N],g[N];
int main(){
	rint n,c,d,q;
	scanf("%d%d%d%d",&n,&c,&d,&q);
	while(q--){
		for(rint i=1;i<=n;i++){
			scanf("%lld",&b[i]);
			b[i]=b[i]*fpow(fpow(i,d),p-2)%p;
		}
		for(rint i=1;i<=n;i++)
			g[i]=fpow(i,p-1+(c-d)%(p-1));
		for(rint i=1;i<=n;i++){
			for(rint j=i+i;j<=n;j+=i){
				g[j]=(g[j]-g[i]+p)%p;
				b[j]=(b[j]-b[i]+p)%p;
			}
			if(!g[i]&&b[i]){
				puts("-1");
				goto E;
			}
			b[i]=b[i]*fpow(g[i],p-2)%p;
		}
		for(rint i=n;i;i--){
			for(rint j=i+i;j<=n;j+=i)
				b[i]=(b[i]-b[j]+p)%p;
		}
		for(rint i=1;i<=n;i++)
			b[i]=b[i]*fpow(fpow(i,d),p-2)%p;
		for(rint i=1;i<=n;i++)
			printf("%lld ",b[i]);
		puts("");E:;
	}
	return 0;
}