UR#5 怎样跑得更快

Posted 2021-05-25 An_Fly

给你整数 \\(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;
}