QOJ 59. Determinant of A+Bz 題解

QOJ # 59. Determinant of A+Bz 題解

起因是有人在聯(lián)考里考了搬了一個需要用到這個的題，本來原題是 \(\mathcal{O}(n^4)\) 就可以過的，但是喪心病狂的出題人把它加強到了 \(\mathcal{O}(n^3)\)。

題意

給定一個 \(n\times n\) 的矩陣，每個位置是一個 \(Ax+B\) 的多項式。你需要求出這個矩陣的行列式，輸出這個 \(n\) 次多項式。對 \(9998244353\) 取模。

\(n\le 500,0\le A_{i,j},B_{i,j} < 998244353\)。

海森堡算法

海森堡可以在 \(\mathcal{O}(n^3)\) 的時間復雜度內(nèi)求解一個矩陣的特征多項式，其中矩陣 \(A\) 的特征多項式為 \(det(I_nx-A)\)。

首先一個暴力做法是用拉格朗日插值求解，即把 \(x\) 帶入 \(0\sim n\) 然后求解行列式，復雜度為 \(\mathcal{O}(n^4)\)。我們發(fā)現(xiàn)復雜度的瓶頸為每次求解行列式，所以可以思考將矩陣變?yōu)橐粋€比較好看的形式，使得每次求行列式可以快速計算。

初等變換

一個矩陣 \(A\) 的初等變換共有三種：

交換第 \(s\) 行和第 \(t\) 行，矩陣表示為將 \(I_n\) 的第 \(s\) 行和第 \(t\) 行互換后的矩陣，記為 \(P_{s,t}\)。\(P_{s,t}A\) 表示交換 \(A\) 的第 \(s\) 行和第 \(t\) 行，\(AP_{s,t}\) 表示交換 \(A\) 的第 \(s\) 列和第 \(t\) 列。
將第 \(s\) 行乘一個數(shù) \(c\)，矩陣表示為將 \(I_n\) 的 \((s,s)\) 位置設(shè)為 \(c\) 的矩陣，記為 \(D_s(c)\)。\(D_s(c)A\) 表示將 \(A\) 的第 \(s\) 行乘 \(c\)，\(AD_s(c)\) 表示將 \(A\) 的第 \(i\) 列乘 \(c\)。
將第 \(s\) 行乘一個數(shù) \(k\) 加到第 \(t\) 行上，矩陣表示為將 \(I_n\) 的 \((t,s)\) 位置設(shè)為 \(k\) 的矩陣，記為 \(T_{s,t}(k)\)。\(T_{s,t}(k)A\) 表示將 \(A\) 的第 \(s\) 行乘 \(k\) 加到第 \(t\) 行上。\(AT_{s,t}(k)\) 表示將 \(A\) 的第 \(t\) 列乘 \(k\) 加到第 \(s\) 列上。

三個矩陣都是可逆矩陣，有 \(P^{-1}_{s,t} = P_{s,t},D^{-1}_s(c) = D_s(c^{-1}),T^{-1}_{s,t}(k) = T_{s,t}(-k)\)。

相似變換

對于 \(n\times n\) 的矩陣 \(A\) 和 \(B\)，當存在可逆矩陣 \(P\) 滿足：

\[B = P^{-1}AP \]

則稱為 \(A\) 和 \(B\) 相似，記變換 \(A\leftrightarrow P^{-1}AP\) 為相似變換。

命題：相似變換不改變矩陣的行列式和特征多項式，即 \(A\) 和 \(P^{-1}AP\) 有相同的行列式和特征多項式。

證明：首先 \(A\) 和 \(P^{-1}AP\) 的行列式相同是好證的，考慮證明特征多項式相同：

\[\begin{aligned} det(I_nx-P^{-1}AP) &= det(I_nxP^{-1}P-P^{-1}AP) \\ &= det(P^{-1}I_nxP-P^{-1}AP) \\ &= det(P^{-1}(I_nx-A)P) \\ &= det(I_nx-A)det(P^{-1})det(P) \\ &= det(I_nx-A) \end{aligned} \]

把 \(P\) 寫成初等矩陣的積，設(shè) \(P = T_1T_2\ldots T_n\)，那么：

\[P^{-1}AP = T_n^{-1}\ldots T_1^{-1}AT_1\ldots T_n \]

于是相似變換相當于每次對 \(A\) 做一個初等行變換，再用這個初等矩陣的逆矩陣對 \(A\) 做一個列變換。三個初等變換分別為：

先交換 \(A\) 的第 \(s\) 行和第 \(t\) 行，再交換 \(A\) 的第 \(s\) 列和第 \(t\) 列。
先將 \(A\) 的第 \(s\) 行乘上 \(c\)，再將 \(A\) 的第 \(t\) 列除以 \(c\)。
先將 \(A\) 的第 \(s\) 行乘 \(k\) 加到第 \(t\) 行上，再將 \(A\) 的第 \(t\) 列乘 \(-k\) 加到第 \(s\) 列上。

上海森堡矩陣

~~跟上海沒有關(guān)系~~

因為相似變換不改變特征多項式，所以我們可以嘗試用對 \(A\) 做相似變換把 \(A\) 消成比較好看的形式。

但是很遺憾，簡單的高斯消元無法通過相似變換將 \(A\) 消成上三角矩陣，因為當你用 \((i,i)\) 去消第 \(j\) 行時，第 \(j\) 列又會對第 \(i\) 列有貢獻，所以不太好做。于是考慮可以將 \(A\) 消成上海森堡矩陣。

上海森堡矩陣的定義為，除主對角線及以上和主對角線下面一格的位置，其余位置都為 \(0\)，即所有滿足 \(i > j+1\) 的 \(A_{i,j} = 0\)。這樣在做消元的時候就是用 \((i+1,i)\) 去消第 \(j\) 行，那么逆矩陣就是用第 \(j\) 列對第 \(i+1\) 列有貢獻，而第 \(i+1\) 列在下一步才會枚舉到，沒有影響。于是就可以在 \(\mathcal{O}(n^3)\) 的復雜度內(nèi)將一個矩陣 \(A\) 通過相似變換消成上海森堡矩陣。

將 \(A\) 消成上海森堡矩陣后，現(xiàn)在的問題就是如何求一個上海森堡矩陣的特征多項式。還是用拉差來計算，每次就是要快速計算一個上海森堡矩陣的行列式（下面這一部分和網(wǎng)上不太相同，沒有用一些比較高深的東西）。

我們考慮普通求行列式的算法，每次用 \((i,i)\) 去消第 \(i+1\sim n\) 行第 \(i\) 列的值。而因為這是上海森堡矩陣，所以只有 \(A_{i+1,i}\) 可能有值，其余 \(j>i+1\) 的 \(j\) 都滿足 \(A_{j,i} = 0\)，所以每次消的時候只用將第 \(i+1\) 行消掉即可，復雜度為 \(\mathcal{O}(n^2)\)。

總的時間復雜度為 \(\mathcal{O}(n^3)\)，于是我們就將特征多項式的求法由 \(\mathcal{O}(n^4)\) 優(yōu)化為了 \(\mathcal{O}(n^3)\)。

這一部分的代碼：

ll det(int n)
{
    ll ans = 1;
    for(int i = 1;i <= n;i++)
    {
        if(!f[i][i])
            if(i+1 <= n&&f[i+1][i])swap(f[i],f[i+1]),ans = -ans;
            else return 0;
        (ans *= f[i][i]) %= mod;
        ll s = f[i+1][i]*inv(f[i][i])%mod;
        for(int j = i;j <= n;j++)
            (f[i+1][j] -= f[i][j]*s) %= mod;
    }
    return (ans+mod)%mod;
}
inline int calc(int n,ll x)
{
    for(int i = 1;i <= n;i++)for(int j = 1;j <= n;j++)
        f[i][j] = ((i==j)*x+b[i][j])%mod;
    return det(n);
}
void Hessenberg(int n)
{
    for(int i = 2;i <= n;i++)
    {
        if(!b[i][i-1])
        {
            int p = 0;
            for(int j = i;j <= n;j++)
                if(b[j][i-1]){p = j;break;}
            if(!p)continue;
            for(int k = 1;k <= n;k++)swap(b[i][k],b[p][k]);
            for(int k = 1;k <= n;k++)swap(b[k][i],b[k][p]);
        }
        for(int j = i+1;j <= n;j++)
        {
            ll s = b[j][i-1]*inv(b[i][i-1])%mod;
            for(int k = 1;k <= n;k++)(b[j][k] -= b[i][k]*s) %= mod;
            for(int k = 1;k <= n;k++)(b[k][i] += b[k][j]*s) %= mod;
        }
    }
    for(int i = 0;i <= n;i++)
    {
        ll s = 1;
        for(int j = 0;j <= n;j++)g[j] = !j;
        for(int j = 0;j <= n;j++)if(i != j)
        {
            s = s*(i-j+mod)%mod;
            for(int k = n;~k;k--)
                g[k] = ((k?g[k-1]:0)+g[k]*(mod-j))%mod;
        }
        s = inv(s)*calc(n,i)%mod;
        for(int j = 0;j <= n;j++)(F[j] += g[j]*s) %= mod;
    }
}

題解

如果我們能把 \(det(Ax+B)\) 的形式轉(zhuǎn)化成 \(det(I_nx+C)\) 的形式，就可以用上面的海森堡算法了。（注意 \(C\) 不取反也行）

首先第一想法是乘一個 \(A^{-1}\)，變成 \(\frac{det(A^{-1}(Ax+B))}{det(A^{-1})} = \frac{det(I_nx+A^{-1}B)}{det(A^{-1})}\)，但是如果 \(A\) 不是滿秩的就無法使用這個做法。

我們考慮直接將 \(A,B\) 兩個矩陣放在一起，然后嘗試將 \(A\) 消成 \(I_n\)。假設(shè)到了第 \(i\) 行，先用前 \(i-1\) 列將第 \(i\) 列的前 \(i-1\) 行都消成 \(0\)，然后找這一列的主元。如果此時 \(A\) 的第 \(i\) 列全是 \(0\)，那么就會有問題。而因為 \(A\) 的第 \(i\) 列全是 \(0\)，所以原矩陣中這一列的元素都不含 \(x\)，我們可以考慮給這一列的所有元素都乘上 \(x\)，最后再將行列式除以 \(x\)。

將一列乘上 \(x\) 相當于將 \(A\) 的第 \(i\) 列賦值成 \(B\) 的第 \(i\) 列然后將 \(B\) 的第 \(i\) 列清空，此時再用前 \(i-1\) 列將第 \(i\) 列的前 \(i-1\) 行都消掉，如果第 \(i\) 列有非零元素就可以找到主元了，而如果這一列仍然全部為 \(0\)，我們就再進行上述操作（消 \(A\) 第 \(i\) 列前 \(i-1\) 行時 \(B\) 也會跟著操作，所以可能會產(chǎn)生新的非零元素）。記錄一個變量表示最后要除多少次 \(x\)，如果中途發(fā)現(xiàn)除的次數(shù) \(>n\) 了，那么最終答案就是 \(0\)。

這樣子做下去，要么是發(fā)現(xiàn)答案為 \(0\)，要么是將 \(A\) 消成了單位矩陣，此時再做海森堡算法即可。整個復雜度依然為 \(\mathcal{O}(n^3)\)。

代碼

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
#define ll long long
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
using namespace std;
const int N = 505,mod = 998244353;
int n,c;
ll a[N][N],b[N][N],f[N][N],g[N],F[N],sum;
inline ll inv(ll x){return x<0?inv((x+mod)%mod):x==1?1:(mod-mod/x)*inv(mod%x)%mod;}
ll det(int n)
{
    ll ans = 1;
    for(int i = 1;i <= n;i++)
    {
        if(!f[i][i])
            if(i+1 <= n&&f[i+1][i])swap(f[i],f[i+1]),ans = -ans;
            else return 0;
        (ans *= f[i][i]) %= mod;
        ll s = f[i+1][i]*inv(f[i][i])%mod;
        for(int j = i;j <= n;j++)
            (f[i+1][j] -= f[i][j]*s) %= mod;
    }
    return (ans+mod)%mod;
}
inline int calc(int n,ll x)
{
    for(int i = 1;i <= n;i++)for(int j = 1;j <= n;j++)
        f[i][j] = ((i==j)*x+b[i][j])%mod;
    return det(n);
}
void Hessenberg(int n)
{
    for(int i = 2;i <= n;i++)
    {
        if(!b[i][i-1])
        {
            int p = 0;
            for(int j = i;j <= n;j++)
                if(b[j][i-1]){p = j;break;}
            if(!p)continue;
            for(int k = 1;k <= n;k++)swap(b[i][k],b[p][k]);
            for(int k = 1;k <= n;k++)swap(b[k][i],b[k][p]);
        }
        for(int j = i+1;j <= n;j++)
        {
            ll s = b[j][i-1]*inv(b[i][i-1])%mod;
            for(int k = 1;k <= n;k++)(b[j][k] -= b[i][k]*s) %= mod;
            for(int k = 1;k <= n;k++)(b[k][i] += b[k][j]*s) %= mod;
        }
    }
    for(int i = 0;i <= n;i++)
    {
        ll s = 1;
        for(int j = 0;j <= n;j++)g[j] = !j;
        for(int j = 0;j <= n;j++)if(i != j)
        {
            s = s*(i-j+mod)%mod;
            for(int k = n;~k;k--)
                g[k] = ((k?g[k-1]:0)+g[k]*(mod-j))%mod;
        }
        s = inv(s)*calc(n,i)%mod;
        for(int j = 0;j <= n;j++)(F[j] += g[j]*s) %= mod;
    }
}
char buf[1<<21],*p1,*p2;
inline int rd()
{
    char c;int f = 1;
    while(!isdigit(c = getchar()))if(c=='-')f = -1;
    int x = c-'0';
    while(isdigit(c = getchar()))x = x*10+(c^48);
    return x*f;
}
int main()
{
    // freopen(".in","r",stdin);
    // freopen(".out","w",stdout);
    n = rd();ll sum = 1;
    for(int i = 1;i <= n;i++)for(int j = 1;j <= n;j++)b[i][j] = rd();
    for(int i = 1;i <= n;i++)for(int j = 1;j <= n;j++)a[i][j] = rd();
    for(int i = 1;i <= n;i++)
    {
        int p = 0;
        while(!p)
        {
            for(int j = i;j <= n;j++)
                if(a[j][i]){p = j;break;}
            if(p)break;
            for(int j = 1;j < i;j++)if(a[j][i])
            {
                ll s = a[j][i]*inv(a[j][j])%mod;
                for(int k = 1;k <= n;k++)
                    (a[k][i] -= a[k][j]*s) %= mod,(b[k][i] -= b[k][j]*s) %= mod;
            }
            c++;for(int j = 1;j <= n;j++)a[j][i] = b[j][i],b[j][i] = 0;
            if(c > n){for(int i = 0;i <= n;i++)printf("%d%c",0," \n"[i==n]);return 0;}
        }
        swap(a[i],a[p]);swap(b[i],b[p]);
        if(i != p)sum = -sum;
        ll s = inv(a[i][i]);(sum *= a[i][i]) %= mod;
        for(int j = 1;j <= n;j++)
            (a[i][j] *= s) %= mod,(b[i][j] *= s) %= mod;
        for(int j = 1;j <= n;j++)if(i != j)
        {
            ll s = a[j][i]*inv(a[i][i])%mod;
            for(int k = 1;k <= n;k++)
                (a[j][k] -= a[i][k]*s) %= mod,(b[j][k] -= b[i][k]*s) %= mod;
        }
    }
    Hessenberg(n);
    for(int i = 0;i <= n;i++)F[i] = i<=n-c?F[i+c]*(sum+mod)%mod:0;
    for(int i = 0;i <= n;i++)printf("%lld%c",F[i]," \n"[i==n]);
    return 0;
}

一些例題

Ascending Matrix - Problem - QOJ.ac：LGV 引理+生成函數(shù)，這題本來 \(\mathcal{O}(n^4)\) 能過，但是放到聯(lián)考里被加強了。
[ABC323G] Inversion of Tree：矩陣樹定理+生成函數(shù)。

posted @ 2025-06-12 11:09 max0810 閱讀(108) 評論(1) 收藏舉報

刷新頁面返回頂部

max0810