1. 程式人生 > >求矩陣的特徵多項式(c++)

求矩陣的特徵多項式(c++)

一個比較慢的做法

  首先你要知道矩陣的特徵多項式是什麼。

  直接消元就可以了。

  時間複雜度:O(n5)O(n5)或O(n4)O(n4)。

一個稍微快一點的做法

  觀察到特徵多項式的次數是nn。

  我們就可以插值。

  具體來說,先求出當x=0…nx=0…n時特徵多項式對應的點值,然後直接用拉格朗日插值插出來。

  時間複雜度:O(n4)O(n4)

一個更快的做法

  有一個性質:相似矩陣的特徵多項式相同。

  所以可以把這個矩陣相似到一個可以快速求特徵多項式的矩陣再求。

  怎麼相似呢?

  AA和BB相似意味著A=PBP−1A=PBP−1

  假設當前矩陣是AA。在消元過程中會左乘一個初等矩陣PP,這時在AA的右邊乘上這個初等矩陣的逆矩陣P−1P−1,得到PAP−1PAP−1。

  顯然這個矩陣和矩陣AA相似。

  但是我們只能把所有i≥j+2i≥j+2的部分消成00,也就是說會消成一個上海森堡矩陣。

  接下來就要快速求一個上海森堡矩陣的特徵多項式。

  上海森堡矩陣長這樣:

⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜a1,1a2,10⋮0a1,2a2,2a3,2⋮0a1,3a2,3a3,3⋮0⋯⋯⋯⋱⋯a1,na2,na3,n⋮an,n⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟(a1,1a1,2a1,3⋯a1,na2,1a2,2a2,3⋯a2,n0a3,2a3,3⋯a3,n⋮⋮⋮⋱⋮000⋯an,n)

  設右下角i×ii×i的矩陣的特徵多項式為fifi。

  考慮對第一列展開。

  那麼肯定是連續的消掉若干個第二行後再消掉第一行。

  即

fi=ai,ifi+1−ai+1,iai,i+1fi+2+ai+1,iai+2,i+1ai,i+2fi+3+⋯fi=ai,ifi+1−ai+1,iai,i+1fi+2+ai+1,iai+2,i+1ai,i+2fi+3+⋯

  注意到只有主對角線上的元素一次的,其他都是常數的,所以前面的係數的次數最多是11。

  這樣可以從後往前求出所有fifi。時間複雜度是O(n3)O(n3)。

  這樣我們就得到了一個時間複雜度是O(n3)O(n3)的優秀做法。

  這個做法有什麼用?

  可以去卡別人。比如用O(n3+n2logm)O(n3+n2log⁡m)的做法卡O(n3logm)O(n3log⁡m)的矩陣快速冪。

程式碼

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<utility>
using namespace std;
typedef long long ll;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
const ll p=998244353;
struct poly
{
    ll a[510];
    int n;
    poly()
    {
        memset(a,0,sizeof a);
        n=0;
    }
    ll &operator [](int x)
    {
        return a[x];
    }
};
int n;
ll a[510][510];
void add(poly &a,poly &b,pll c)
{
    a.n=max(a.n,b.n+bool(c.first));
    int i;
    for(i=0;i<=a.n;i++)
    {
        a[i]=(a[i]+b[i]*c.second)%p;
        if(i)
            a[i]=(a[i]+b[i-1]*c.first)%p;
    }
}
ll fp(ll a,ll b)
{
    ll s=1;
    for(;b;b>>=1,a=a*a%p)
        if(b&1)
            s=s*a%p;
    return s;
}
void gao1()
{
    int i,j,k;
    for(i=1;i<=n;i++)
    {
        for(j=i+1;j<=n;j++)
            if(a[i][j])
                break;
        if(j>n)
            continue;
        if(j!=i+1)
        {
            for(k=i;k<=n;k++)
                swap(a[i+1][k],a[j][k]);
            for(k=1;k<=n;k++)
                swap(a[k][i+1],a[k][j]);
        }
        ll e=fp(a[i+1][i],p-2);
        for(j=i+2;j<=n;j++)
            if(a[j][i])
            {
                ll v=e*a[j][i]%p;
                for(k=i;k<=n;k++)
                    a[j][k]=(a[j][k]-a[i+1][k]*v)%p;
                for(k=1;k<=n;k++)
                    a[k][i+1]=(a[k][i+1]+a[k][j]*v)%p;
            }
    }
}
pll c[510][510];
poly f[510];
pll operator *(pll a,pll b)
{
    return pll((a.first*b.second+a.second*b.first)%p,a.second*b.second%p);
}
pll operator *(pll a,ll b)
{
    return pll(a.first*b%p,a.second*b%p);
}
void gao2()
{
    f[n+1][0]=1;
    int i,j;
    for(i=n;i>=1;i--)
    {
        pll v(0,1);
        for(j=i+1;j<=n+1;j++)
        {
            add(f[i],f[j],v*c[i][j-1]*((j-i)&1?1:-1));
            v=v*c[j][j-1];
        }
    }
}
int main()
{
    scanf("%d",&n);
    int i,j;
    for(i=1;i<=n;i++)
        for(j=1;j<=n;j++)
            scanf("%lld",&a[i][j]);
    gao1();
    for(i=1;i<=n;i++)
        for(j=1;j<=n;j++)
        {
            c[i][j].second=-a[i][j];
            if(i==j)
                c[i][j].first=1;
        }
    gao2();
    for(i=0;i<=n;i++)
        printf("%lld ",(f[1][i]+p)%p);
    printf("\n");
    return 0;
}