Jordan形$O(n^4)$算法
Computing Jordan canonical form
There are better algorithms for computing Jordan form. For example, An \(O(n^3)\) Algorithm for the Frobenius Normal Form Storjohann 1998 is an easy faster way.
But it is hard to compute transfer matrix using these faster methods.
Here we give the proof and algorithm simultaneously
Definition 1 root subspace for \(\lambda\) is the subspace satisfying \(\lambda\) \(\exist k, (\mathcal{A}-\lambda\mathcal{I})^k\alpha=0\)
Lemma 1 The sum of different root subspaces is direst sum
proof
\(f(x)=(x-\lambda_1)^t\) is the characteristic polynomial's factor of \(\lambda_1\) , \(g(x)\) is the remaining part
Now \(f(x)g(x)\) is annihilator for \(W_{\lambda _1}\) (Hamilton-Cayley)
\(g(x)\) cannot annihilate any nonzero vector in \(W_{\lambda _1}\) , because if \(g\) can , \(\gcd(f,g)\) can by bezout's formula, but \(\gcd(f(x),g(x))=1\)
so \(f(x)\) is annihilator.
to prove \(W_{\lambda_1}\oplus W_{\lambda_2}\dots\oplus W_{\lambda_t}\) we only need to prove vector \(0\) have only one representation.
say \(0=u_1+u_2+\dots+u_t\) . applying \(g(\mathcal{A})\) to both side yield \(g(\mathcal{A})u_1=0\) .
and because \(\gcd(f(x),g(x))=1\) and \(f(\mathcal{A})u_1=0\) , \(\mathcal{I}u_1=u_1=0\)
\(\square\)
Definition 2 cyclic subspace \(C(\alpha)\) is the subspace generated by \(\alpha \ ,\mathcal{A}\alpha,\ \mathcal{A}^2\alpha \dots\)
it's easy to prove \(\alpha \ ,\mathcal{A}\alpha,\ \mathcal{A}^2\alpha \dots\) are independent. but it is not needed in the sequel
Lemma 2 \(V\) is the direct sum of some cyclic subspace for nilpotent linear transformation \(\mathcal{A}\)
proof
we give a constructive and algorithmic proof.
Let the degree of minimal polynomial be \(d\). (The minmimum \(d\) satisfying \(A^d=0\) )
we claim that \(V=\oplus_{k=1}^d ker(A^k)/ker(A^{k-1})\) (\(\star\)) because
Here when we talk about a vector in \(ker(A^k)/ker(A^{k-1})\) we mean a representative in the whole space \(V\).
The Algorithm is as follows.
Denote \(X\) the set containing a basis of \(V\) .
start with \(X\) empty.
Insert the basis of \(ker(A^d)/ker(A^{d-1})\) into \(X\) , denote them \(\alpha_1,\alpha_2,\dots\)
take out the elements in \(X\) now (they are \(\alpha_1,\alpha_2,\dots\) ), and apply \(A\) on them.
Now \(A\alpha_1,A\alpha_2,\dots\) are vectors in \(ker(A^{d-1})/ker(A^{d-2})\) . and they(their projection image into quotient space) are linearly independent in \(ker(A^{d-1})/ker(A^{d-2})\) because if
then
that is
\(\{\alpha_i\}\) is a basis so all \(\mu_i=0\)
we can enlarge \(A\alpha_1,A\alpha_2,\dots\) to a basis of \(ker(A^{d-1})/ker(A^{d-2})\), denote the new vectors \(\beta_1,\beta_2,\dots\)
Insert \(A\alpha_1,A\alpha_2,\dots \beta_1,\beta_2,\dots\) into \(X\)
because of (\(\star\)) \(\alpha_1,\alpha_2,\dots\) are independent with newly inserted vectors in \(V\)
now \(X\) is \(\{\alpha_1,\alpha_2,\dots A\alpha_1,A\alpha_2,\dots \beta_1,\beta_2,\dots\}\)
take out the elements in \(X\) inserted in the previous step (they are \({A\alpha_1,A\alpha_2,\dots \beta_1,\beta_2,\dots}\)
), and apply \(A\) on them.
the result (\({A^2\alpha_1,A^2\alpha_2,\dots A\beta_1,A\beta_2,\dots}\)) is linearly independent in \(ker(A^{d-2})/ker(A^{d-3})\) so we can enlarge them to a basis.
continue this process.
Ultimately, \(X\) is a basis of \(V\) because of (\(\star\))
\(\square\)
Note that to computationally determine whether some vectors \(\alpha_1,\alpha_2,\dots\) are linearly independent in \(ker(A^k)/ker(A^{k-1})\) we need to just determine whether \(A^{k-1}\alpha_1,A^{k-1}\alpha_2,\dots\) are linearly independent in \(V\) . You can demonstrate this by write down the definition. (again \(\alpha_i\) are representatives in V, not vectors in qutient space)
Theorem There is a basis under which the matrix is in Jordan form.
proof
For every eigenvalue, run algorithm of lemma 2. And the union of vectors is a Jordan basis.
Furthermore, it is easy to see that the number of different size of Jordan square is determined by the solution of \((A-\lambda^kI)X=0\) so it is determined by \(rank(A-\lambda I)^k\) . Which is invariant for similar relation.
So Jordan form is canonical form of similar matrices. \(\square\)
Corollary Jordan form and its transfer matrix can be computed in \(O(n^4)\).
proof
Computing the solution for \(A^kX=0\) need \(Cn^4\) operations if implemented using ordinary method.
Computing the chain of a vector \(\alpha\) (the space \(C(\alpha)\)) need \(C\sum d_i\times n^2=Cn^3\). Maintaining the kernel (because we need to determine the linear dependency in the quotient space) need \(C\sum d_i^3=O(n^3)\)
So the complexity is \(O(n^4)\) \(\square\)
The code below needs eigenvalues.
#include<bits/stdc++.h>
#define de fprintf(stderr,"on-%d\n",__LINE__)
#define pb push_back
#define fr(i,a,b) for(int i(a),end##i(b);i<=end##i;i++)
#define fd(i,a,b) for(int i(a),end##i(b);i>=end##i;i--)
#define reduce Reduce
#define apply Apply
using namespace std;
const int N=30;
int m;
struct frac{
int q,p;
frac(){
q=0;p=1;
}
void out(){
if(p==1)printf("%7d ",q);
else printf("%3d/%3d ",q,p);
}
};
frac zero;
frac make(int x,int y){
frac A;
A.q=x;A.p=y;
return A;
}
frac r[N];
int d[N];
frac V(int x){
return make(x,1);
}
frac read(){
int x;scanf("%d",&x);
return V(x);
}
int gcd(int x,int y){
if(!x)return y;
return gcd(y%x,x);
}
frac reduce(frac &A){
int g=gcd(A.p,A.q);
A.p/=g;A.q/=g;
if(A.p<0)A.p=-A.p,A.q=-A.q;
return A;
}
frac reduce(frac A){
int g=gcd(A.p,A.q);
A.p/=g;A.q/=g;
if(A.p<0)A.p=-A.p,A.q=-A.q;
return A;
}
frac operator +(const frac &A,const frac &B){
return reduce(make(A.q*B.p+B.q*A.p,A.p*B.p));
}
frac operator -(const frac &A,const frac &B){
return reduce(make(A.q*B.p-B.q*A.p,A.p*B.p));
}
frac operator *(const frac &A,const frac &B){
assert(A.p);
return reduce(make(A.q*B.q,A.p*B.p));
}
frac operator /(const frac &A,const frac &B){
return reduce(make(A.q*B.p,A.p*B.q));
}
struct matrix{
frac A[N][N];
frac* operator [](int x){return A[x];}
void out(){
fr(i,1,m){
fr(j,1,m)A[i][j].out();
printf("\n");
}
}
};
matrix I;
void init(){
fr(i,1,m)fr(j,1,m)I[i][j]=V(i==j);
}
matrix operator *(matrix A,matrix B){
matrix C;
fr(i,1,m)fr(j,1,m){
C[i][j]=V(0);
fr(k,1,m)C[i][j]=C[i][j]+A[i][k]*B[k][j];
}
return C;
}
matrix operator *(frac v,matrix A){
fr(i,1,m)fr(j,1,m)A[i][j]=v*A[i][j];
return A;
}
matrix operator +(matrix A,matrix B){
fr(i,1,m)fr(j,1,m)A[i][j]=A[i][j]+B[i][j];
return A;
}
matrix operator -(matrix A,matrix B){
fr(i,1,m)fr(j,1,m)A[i][j]=A[i][j]-B[i][j];
return A;
}
vector<frac> apply(matrix &A,const vector<frac>& B){
vector<frac> C;C.resize(m+1);
fr(i,1,m)fr(k,1,m){
C[i]=C[i]+A[i][k]*B[k];
}
return C;
}
bool p[N];
int to[N];
void rowreduce(matrix &A,int n,int m,int op=0){
fr(i,1,m)to[i]=p[i]=0;
fr(i,1,n){
int s=0;
fr(j,1,m)if(!p[j]&&A[i][j].q){
s=j;break;
}
if(!s)continue;
to[i]=s;
p[s]=1;
fr(j,op*i+1,n)if(j!=i&&A[j][s].q){
frac t=A[j][s]/A[i][s];
fr(k,1,m)A[j][k]=A[j][k]-t*A[i][k];
}
}
}
void swap(frac* A,frac *B){
fr(i,1,m)swap(A[i],B[i]);
}
matrix inverse(matrix A){
matrix B=I;
fr(i,1,m){
int s=0;
fr(j,i,m)if(A[j][i].q){
s=j;break;
}
assert(s);
swap(A[i],A[s]);swap(B[i],B[s]);
const frac t=A[i][i];
fr(j,1,m){
A[i][j]=A[i][j]/t;
B[i][j]=B[i][j]/t;
}
fr(j,1,m)if(j!=i&&A[j][i].q){
const frac t=A[j][i];
fr(k,1,m){
A[j][k]=A[j][k]-t*A[i][k];
B[j][k]=B[j][k]-t*B[i][k];
}
}
}
return B;
}
vector<vector<frac> > solve(matrix A){
vector<vector<frac> >W;
rowreduce(A,m,m,0);
fr(i,1,m)if(!p[i]){
vector<frac> T;T.resize(m+1);
T[i]=V(1);
fr(j,1,m)if(to[j]&&A[j][i].q){
T[to[j]]=zero-A[j][i]/A[j][to[j]];
}
W.push_back(T);
}
return W;
}
bool noempty(frac* A){
fr(i,1,m)if(A[i].q)return 1;
return 0;
}
int main(){
#ifdef pig
freopen("pig.in","r",stdin);
freopen("pig.out","w",stdout);
#endif
cin>>m;
init();
matrix A;
fr(i,1,m)fr(j,1,m)A[i][j]=read();
int s;cin>>s;
fr(i,1,s){
r[i]=read();
scanf("%d",&d[i]);
}
matrix Q;int top=0;
fr(i,1,s){
matrix B=A-r[i]*I,C=B;
vector<vector<vector<frac> > > V;
vector<vector<frac> > P;
fr(j,1,d[i]){
vector<vector<frac> > W=solve(C);
matrix A;
fr(k,0,P.size()-1){
fr(l,1,m)A[k+1][l]=P[k][l];
}
fr(k,0,W.size()-1){
fr(l,1,m){
A[k+P.size()+1][l]=W[k][l];
}
}
rowreduce(A,P.size()+W.size(),m,1);
vector<vector<frac> >T;
fr(k,P.size()+1,W.size()+P.size())if(noempty(A[k])){
vector<frac> TT;TT.resize(m+1);
fr(l,1,m)TT[l]=A[k][l];
T.pb(TT);
}
V.pb(T);
P=W;
C=C*B;
}
vector<vector<frac> >ans;
matrix now;int tot=0;
fd(j,V.size()-1,0){
fr(k,0,V[j].size()-1){
vector<frac> T; T.resize(m+1);
fr(l,1,m){
T[l]=V[j][k][l];
}
fr(l,0,j-1)T=apply(B,T);
tot++;
fr(l,1,m)now[tot][l]=T[l];
}
rowreduce(now,tot,m,1);
fr(k,tot-V[j].size()+1,tot)if(noempty(now[k])){
vector<frac> T;T.resize(m+1);
fr(l,1,m)T[l]=V[j][k-(tot-V[j].size()+1)][l];
fr(t,0,j){
ans.pb(T);
T=apply(B,T);
}
}
}
fr(j,0,ans.size()-1){
top++;
fr(k,1,m)Q[k][top]=ans[j][k];
}
}
assert(top==m);
Q.out();
printf("\n");
(inverse(Q)*A*Q).out();
return 0;
}
