求矩陣的特徵多項式(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+n2logm)的做法卡O(n3logm)O(n3logm)的矩陣快速冪。
程式碼
#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;
}