题目来源: 2022 ICPC 华为上海训练营 Day3: G. Matrices and Determinants

题面

一共 \(T\) 组数据.

给定 \(n\times n\) 矩阵 \(\newcommand{\mathbs}[1]{\boldsymbol{ #1 }}\newcommand{\mbs}{\mathbs}\newcommand{\mbb}{\mathbb}\mbs A\in M_n(\mbb Z)\), 其所有元素均为整数. 请求出矩阵 \(\mbs B,\mbs C\in M_n(\mbb Z)\), 所有元素均为整数, 满足: \(\mbs B\mbs C=\mbs A\)\(\newcommand{\abs}[1]{\left| #1 \right|}\abs{\mbs B}=\abs{\mbs C}\ne0\). 如无解, 输出 No. 如有多组满足题意的解, 请输出任意一种.

数据范围

\(1\le T\le 10000\).

\(1\le n\le 4\), \(\mbs A\) 中的所有元素均为整数,且元素 \(A_{ij}\) 满足 \(-10\le A_{ij}\le10\).

题解

由于 \(\mbs {BC}=\mbs A\), 故 \(\abs{\mbs A}=\abs{\mbs{B}}\abs{\mbs C}=\abs{\mbs B}^2\). 由于 \(\mbs {A},\mbs B,\mbs C\) 的元素均为整数, 其行列式也为整数. 故有解的必要条件为: \(\abs{\mbs A}\) 是完全平方数. 且由于 \(\abs{\mbs B}\ne0\), 可知 \(\abs{\mbs A}\ne0\), 故有解的另一个必要条件为: \(\mbs A\) 是非奇异矩阵.

先求 \(\mbs A\) 的行列式. 做 Gauss 消元.

注意到以下事实: 在进行一项初等行变换时, 设变换前的矩阵为 \(\mbs X\), 变换后的矩阵为 \(\mbs Y\), 则有 \[ \mbs Y=\mbs T\mbs X \] 其中 \(\mbs T\) 是对单位阵 \(\mbs I\) 做相同的初等行变换得到的结果. 对于三种初等行变换, \(\mbs T\) 的逆矩阵有三种情况:

  1. 将第 \(k\) 行乘 \(s\), 变换矩阵 \(\mbs T=:\mbs M(k,s)\) 的逆矩阵是 \(\mbs M(k,s^{-1})\).
  2. 将第 \(k\) 行与第 \(l\) 行交换, 变换矩阵 \(\mbs T=:\mbs E(k,l)\) 的逆矩阵是 \(\mbs E(l,k)=\mbs T\), 即逆矩阵为原矩阵.
  3. 将第 \(k\) 行的 \(s\) 倍加到第 \(l\) 行, 变换矩阵 \(\mbs T=:\mbs P(k,s,l)\) 的逆矩阵是 \(\mbs P(k,-s,l)\).

此时有 \(\mbs X=\mbs T^{-1}\mbs Y\).

若在进行初等行变换时, 进行如下限制, 在 \(\mbs X\) 是整数矩阵的情况下, 可以确保 \(\mbs T,\mbs T^{-1}\)\(\mbs Y\) 不会产生非整数:

  1. 进行第 1 种初等行变换时, \(s\) 只取 \(-1\).
  2. 第 2 种初等行变换没有限制.
  3. 进行第 3 种初等行变换时, \(s\) 只取整数.

考虑 Gauss 消元, 在上述限制条件下将 \(\mbs A\) 转化为上三角矩阵的过程中, 需使用辗转相除法. 将 \(\mbs A\) 的第 \(j\) 行的第 \(i\) 列化为 \(0\), 若直接取变换矩阵为 \(\mbs P(i,-\frac{\mbs A[j;i]}{\mbs A[i;i]},j)\), 则不能保证 \(-\frac{\mbs A[j;i]}{\mbs A[i;i]}\) 是整数, 故取变换矩阵为 \[ \newcommand{\brac}[1]{\left( #1 \right)}\newcommand{\flor}[1]{\left\lfloor #1 \right\rfloor}\mbs P\brac{i,-\flor{\frac{\mbs A[j;i]}{\mbs A[i;i]}},j} \] 交换第 \(i\) 行和第 \(j\) 行后再重复进行上述过程, 直至 \(\mbs A[j;i]=0\) 为止.

\(\mbs A\) 经过 \(k\) 次初等行变换后, 得到上三角矩阵 \(\mbs A_k\), 则有 \[ \mbs A=\mbs T_1^{-1}\mbs T_2^{-1}\cdots\mbs T_k^{-1}\mbs A_k=:\mbs J\mbs A_{k} \] 对于 Gauss 消元, 化为上三角矩阵的过程中, 只用到上述的第 2,3 种初等行变换, 此时 \(\mbs J\) 的行列式为 \(1\)\(-1\). 设进行了 \(k_2\) 次第 2 种初等行变换, 则 \(\abs{\mbs J}=(-1)^{k_2}\). \(\mbs A_k\) 是上三角矩阵, 行列式为其主对角线元素的乘积. 此时可以求得 \(\mbs A\) 的行列式了.

若满足必要条件:

由于上三角矩阵的乘积仍为上三角矩阵, 故考虑将 \(\mbs A_k\) 分解为两个上三角矩阵的乘积. 我们希望: \(\mbs A_k=\mbs D\mbs C\), 其中 \(\abs{\mbs C}=\abs{\mbs J}\abs{\mbs D}\), 这样令 \(\mbs B=\mbs {JD}\) 即满足题意.

设: \[ \newcommand{\pmat}[1]{\begin{pmatrix} #1 \end{pmatrix}}\mbs A_k=\pmat{a_{11}&a_{12}&\cdots&a_{1n}\\&a_{22}&\cdots&a_{2n}\\&&\ddots&\vdots\\&&&a_{nn}} \]

\[ \mbs D=\pmat{b_{11}&b_{12}&\cdots&b_{1n}\\&b_{22}&\cdots&b_{2n}\\&&\ddots&\vdots\\&&&b_{nn}} \]

\[ \mbs C=\pmat{c_{11}&c_{12}&\cdots&c_{1n}\\&c_{22}&\cdots&c_{2n}\\&&\ddots&\vdots\\&&&c_{nn}} \]

根据矩阵乘法的定义和上三角矩阵的特点, 有 \[ a_{ii}=b_{ii}c_{ii}\tag1 \]

\(i\ne j\) 时有 \[ a_{ij}=\sum_{k=i}^jb_{ik}c_{kj}=b_{ii}c_{ij}+\sum_{k=i+1}^{j-1}b_{ik}c_{kj}+b_{ij}c_{jj}\tag2 \]

\(\newcommand{\lrac}[1]{\left\{ #1 \right\}}\mbs b_i:=\lrac{b_{1,1+i},b_{2,2+i},\cdots,b_{n-i,n}}\)\(b_{1,1+i}\) 所在的对角线的元素, 同理记 \(\mbs c_i:=\lrac{c_{1,1+i},c_{2,2+i},\cdots,c_{n-i,n}}\). 则可知当 \(\mbs b_0,\mbs b_1,\cdots,\mbs b_{j-i-1}\)\(\mbs c_0,\mbs c_1,\cdots,\mbs c_{j-i-1}\) 均确定后, \[ \sum_{k=i+1}^{j-1}b_{ik}c_{kj} \] 是可以求出的常数, 记做 \(K\). 此时 \((2)\) 式变为 \[ c_{jj}b_{ij}+b_{ii}c_{ij}=a_{ij}-K\tag3 \] 其中 \(b_{ii},c_{jj},a_{ij},K\) 均为已知量, 而 \(b_{ij},c_{ij}\) 是未知量. 此时若满足 \[ \gcd (b_{ii},c_{jj})\mid a_{ij}-K \] 可以根据 \((3)\) 式, 利用根据 Bézout 定理求解 \(b_{ij},c_{ij}\). 此时, 能求出 \(b_{ij},c_{ij}\)充分条件为: \(\gcd(b_{ii},c_{jj})=1\).

根据上述充分条件, 构造 \(\mbs b_0,\mbs c_0\), 整理条件如下:

  1. \(\abs{\mbs C}=\sqrt{\abs{\mbs{A}}}\).
  2. \(a_{ii}=b_{ii}c_{ii}\), \(\forall i\le n\).
  3. \(\gcd(b_{ii},c_{jj})=1\), \(\forall i<j\le n\).

构造的基本思路为: 根据条件 3, \(j\) 越大, \(c_{jj}\) 的限制越多. 故基本思路为, 将尽可能多的 \(\sqrt {\abs{\mbs A}}\) 的因子安排到 \(\mbs c_0\) 中靠前的位置. 对于 \(c_{ii}\), 由于限制条件 1,2, 有不等式 \[ c_{ii}\le\gcd\brac{a_{ii},\frac{\sqrt{\abs{\mbs A}}}{c_{11}c_{22}\cdots c_{i-1,i-1}}} \] 现考虑直接取 \[ c_{ii}=\gcd\brac{a_{ii},\frac{\sqrt{\abs{\mbs A}}}{c_{11}c_{22}\cdots c_{i-1,i-1}}}\tag4 \] 并通过条件 2 构造 \(\mbs b_{0}\), 检验是否满足其他条件. \[ \newcommand{\algn}[1]{\begin{aligned} #1 \end{aligned}}\algn{\abs{\mbs C}=c_{11}c_{22}\cdots c_{nn}&=c_{11}c_{22}\cdots c_{n-1,n-1}\gcd\brac{a_{nn},\frac{\sqrt{\abs{\mbs A}}}{c_{11}c_{22}\cdots c_{n-1,n-1}}}\\&=\gcd({a_{nn}c_{11}c_{22}\cdots c_{n-1,n-1},\textstyle{\sqrt{\abs{\mbs A}}}})\\&=\gcd(a_{nn}\gcd(a_{n-1,n-1}c_{11}c_{22}\cdots c_{n-2,n-2},\textstyle{\sqrt{\abs{\mbs A}}}),\textstyle{\sqrt{\abs{\mbs A}}})\\&=\gcd(a_{nn}a_{n-1,n-1}c_{11}c_{22}\cdots c_{n-2,n-2},a_{nn}\textstyle{\sqrt{\abs{\mbs A}}},\textstyle{\sqrt{\abs{\mbs A}}})\\&=\gcd(a_{nn}a_{n-1,n-1}c_{11}c_{22}\cdots c_{n-2,n-2},\textstyle{\sqrt{\abs{\mbs A}}})\\&=\gcd(a_{nn}\cdots a_{11},\textstyle{\sqrt{\abs{\mbs A}}})=\gcd(\abs{\mbs A},\textstyle{\sqrt{\abs{\mbs A}}})=\textstyle{\sqrt{\abs{\mbs A}}} } \] 上式表明了构造满足条件 1.

\(g_{ij}:=\gcd(b_{ii},c_{jj})\ne 1\), 则由于 \[ g_{ij}\mid c_{jj}\mid\frac{\sqrt{\abs{\mbs A}}}{c_{11}c_{22}\cdots c_{j-1,j-1}}\mid\frac{\sqrt{\abs{\mbs A}}}{c_{11}c_{22}\cdots c_{ii}} \]\(g_{ij}\mid b_{ii}=a_{ii}/c_{ii}\). 此时取 \(c_{ii}'=c_{ii}g_{ij}\), 可知 \[ c_{ii}'=c_{ii}g_{ij}\mid c_{ii}b_{ii}=a_{ii} \]

\[ c_{ii}'=c_{ii}g_{ij}\mid c_{ii}\frac{\sqrt{\abs{\mbs A}}}{c_{11}c_{22}\cdots c_{ii}}=\frac{\sqrt{\abs{\mbs A}}}{c_{11}c_{22}\cdots c_{i-1,i-1}} \]

\(c_{ii}'\)\(a_{ii}\)\(\sqrt{\abs{\mbs A}}/\brac{c_{11}c_{22}\cdots c_{i-1,i-1}}\) 的公因子, 故 \[ c_{ii}<c_{ii}'\le \gcd\brac{a_{ii},\frac{\sqrt{\abs{\mbs A}}}{c_{11}c_{22}\cdots c_{i-1,i-1}}} \]\((4)\) 式矛盾. 故构造满足条件 3.

故按此方法构造满足条件: 构造 \(\mbs b_0,\mbs c_0\) 后, 依次求解 \((3)\) 式构造 \(\mbs b_1,\mbs c_1,\mbs b_2,\mbs c_2,\cdots,\mbs b_{n-1},\mbs c_{n-1}\), 即可构造出 \(\mbs D,\mbs C\). 取 \(\mbs B=\mbs {JD}\) 即满足题意.

复杂度分析

Gauss 消元法需进行 \(O(n^2)\) 步, 每步需要进行 \(O(\log C)\) 次初等行变换, 总复杂度为 \(O(n^3\log C)\).

构造 \(\mbs c_0\) 的复杂度为 \(O(n\log C)\). 求解一个 \(b_{ij},c_{ij}\) 时, 求 \(K\) 的复杂度为 \(O(n)\), Bézout 定理的复杂度为 \(O(\log C)\). 共有 \(O(n^2)\)\(b_{ij},c_{ij}\). 故求解 \(\mbs D,\mbs C\) 的复杂度为 \(O(n^3)\).

总复杂度为 \(O(n^3\log C)\)​.

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
struct matrix
{
ll n;
ll val[4][4];
matrix operator*(const matrix &b) const
{
matrix ans;
ans.n = n;
for (ll i = 0; i < n; ++i)
for (ll j = 0; j < n; ++j)
{
ans.val[i][j] = 0;
for (ll k = 0; k < n; ++k)
ans.val[i][j] += val[i][k] * b.val[k][j];
}
return ans;
}
};
ll gcd(ll a, ll b) { return b ? gcd(b, a % b) : a; }
ll exgcd(ll a, ll b, ll &x, ll &y)
{
if (!b) { x = 1; y = 0; return a; }
ll d = exgcd(b, a % b, y, x);
y = y - a / b * x;
return d;
}
int main()
{
#ifndef pasokon
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
#endif
matrix ZERO;
for (ll i = 0; i < 4; ++i)
for (ll j = 0; j < 4; ++j)
ZERO.val[i][j] = 0;
matrix I = ZERO;
for (ll i = 0; i < 4; ++i)
I.val[i][i] = 1;
ll T;
cin >> T;
while (T--)
{
matrix A, B = I;
cin >> A.n;
ll n = A.n;
B.n = I.n = ZERO.n = n;
for (ll i = 0; i < n; ++i)
for (ll j = 0; j < n; ++j)
cin >> A.val[i][j];
matrix SA = A;
ll sum = 1;
// Gauss 消元
for (ll i = 0, j = 0; i < n; ++i, ++j)
{
ll ni = i;
if (A.val[i][j] == 0)
for (ni = i + 1; ni < n; ++ni)
if (A.val[ni][j])
{
for (ll s = j; s < n; ++s)
swap(A.val[i][s], A.val[ni][s]);
sum *= -1;
matrix S = I;
S.val[i][i] = S.val[ni][ni] = 0;
S.val[i][ni] = S.val[ni][i] = 1;
B = B * S;
break;
}
if (ni >= n)
{
sum = 0;
break;
}
for (ll ni = i + 1; ni < n; ++ni)
{
matrix S = I;
S.val[i][i] = S.val[ni][ni] = 0;
S.val[i][ni] = S.val[ni][i] = 1;
while (A.val[ni][j])
{
ll tms = A.val[i][j] / A.val[ni][j];
for (ll s = j; s < n; ++s)
A.val[i][s] -= tms * A.val[ni][s];
for (ll s = j; s < n; ++s)
swap(A.val[i][s], A.val[ni][s]);
sum *= -1;
matrix D = I;
D.val[ni][i] = tms;
B = B * S * D;
}
}
sum *= A.val[i][j];
}
ll sq = round(sqrt(sum));
if (sq * sq != sum || sum == 0)
{
cout << "No\n";
continue;
}
cout << "Yes\n";
matrix C = ZERO, D = ZERO;
// 安排对角线
for (ll j = n - 1; j >= 0; --j)
{
ll g = gcd(abs(A.val[j][j]), sq);
C.val[j][j] = g;
D.val[j][j] = A.val[j][j] / g;
sq /= g;
}
// exgcd 求其他元素
for (ll d = 1; d < n; ++d)
{
// d = j - i
for (ll i = 0, j = d + i; j < n; ++i, ++j)
{
ll r = A.val[i][j];
for (ll k = i + 1; k < j; ++k)
r -= C.val[i][k] * D.val[k][j];
// C.val[i][i] * D.val[i][j] + C.val[i][j] * D.val[j][j] = r;
ll g = exgcd(C.val[i][i], D.val[j][j], D.val[i][j], C.val[i][j]);
D.val[i][j] *= (r / g), C.val[i][j] *= (r / g);
}
}
B = B * C;
ll check = 1;
for (ll i = 0; i < n; ++i)
check *= D.val[i][i];
if (check == 1)
{
matrix II = I;
II.val[0][0] = -1;
B = B * II;
D = II * D;
}
for (ll i = 0; i < n; ++i)
{
for (ll j = 0; j < n; ++j)
cout << (j ? " " : "") << B.val[i][j];
cout << '\n';
}
for (ll i = 0; i < n; ++i)
{
for (ll j = 0; j < n; ++j)
cout << (j ? " " : "") << D.val[i][j];
cout << '\n';
}
}
return 0;
}

花絮

笔者于训练赛开始后4小时59分钟(最后一分钟)通过本题,且获得了本题的最快解题奖,也是全场唯一的通过。