2022 ICPC 华为上海训练营 Day3G 题解

题目来源: 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)\)​.

代码实现

#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分钟（最后一分钟）通过本题，且获得了本题的最快解题奖，也是全场唯一的通过。