//基本(タテ)ベクトルは次のように入力できる。

e1:=matrix([1,0,0,0]):
e2:=matrix([0,1,0,0]):
e3:=matrix([0,0,1,0]):
e4:=matrix([0,0,0,1]):

// A,B は基本ベクトルを次のように並べたものである。

A:=e2.e3.e4.e1:
B:=e2.e1.e4.e3:
print ("A=",A,"B=",B);

v:=e1+e2+e3+e4:

//ひとまず A,B の基底を変換する。
Q:=v.e2.e3.e4:
A1:=Q^(-1)*A*Q:
B1:=Q^(-1)*B*Q:

// 射影

P1:=matrix([[1,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]):  

// P1 の平均を計算する。和をとるには $記号が重宝する。

P1M:=1/8*_plus(A1^i*B1^j*P1*((A1^i*B1^j)^(-1))$ i=0..3 $ j=0..1):

print("P1=",P1,"P1M=",P1M);

// P1M の像と核の基底の計算をするわけだが、
// 射影の場合にはこれは固有ベクトルを求めても同じである。

RESULT2:=linalg::eigenvectors(P1M):

//RESULT2 は固有値、重複度、固有ベクトルの順に並んでいるので、
//ここから必要なデータを得る。

w1:=RESULT2[1][3][1]:
w2:=RESULT2[1][3][2]:
w3:=RESULT2[1][3][3]:
w0:=RESULT2[2][3][1]:

// Q1 は P1M の像と核との基底を並べたもの。

Q1:=w0.w1.w2.w3:

print(
"(Q1^(-1))*A1*Q1=",
(Q1^(-1))*A1*Q1,
"(Q1^(-1))*B1*Q1=",
(Q1^(-1))*B1*Q1);