% an inclusion to diagonalize a matrix

% all  matrices of form [[a00 a01] [a10 a11]]

/E 2.718281828 def

/PI 3.141592654 def

% 2 x 2 matrix utilities -------------------------------------------

% 2 x 2 determinant

/det {
1 dict begin
/A exch def
A 0 get 0 get
A 1 get 1 get
mul
A 0 get 1 get
A 1 get 0 get
mul
sub 
end
} def

% 2 x 2 inverse

/trace {
1 dict begin
/A exch def
A 0 get 0 get
A 1 get 1 get
add
end
} def

/inverse {
2 dict begin
/A exch def
/d A det def
[
[A 1 get 1 get d div
A 0 get 1 get neg d div]
[A 1 get 0 get neg d div
A 0 get 0 get d div]
]
end
} def

% A c -> Ac

/scalarmul {
pstack
2 dict begin
/c exch def
/A exch def

[

/row A 0 get def

[ row 0 get c mul
row 1 get c mul ]

/row A 1 get def

[ row 0 get c mul
row 1 get c mul ]

]

end
} def

% replaces A B by AB

/matrixmul {
2 dict begin
/B exch def
/A exch def
[
	[
	A 0 get 0 get B 0 get 0 get mul
	A 0 get 1 get B 1 get 0 get mul
	add
	A 0 get 0 get B 0 get 1 get mul
	A 0 get 1 get B 1 get 1 get mul
	add
	]
	[
	A 1 get 0 get B 0 get 0 get mul
	A 1 get 1 get B 1 get 0 get mul
	add
	A 1 get 0 get B 0 get 1 get mul
	A 1 get 1 get B 1 get 1 get mul
	add
	]
]
end
} def

% replaces A v by Av

/matrixvectormul {
2 dict begin
/v exch def
/A exch def
[
A 0 get 0 get v 0 get mul
A 0 get 1 get v 1 get mul
add
A 1 get 0 get v 0 get mul
A 1 get 1 get v 1 get mul
add
]
end
} def

% - type -----------------------------------------------------

% first one fixes type Real, Unipotent, Complex
% input = a matrix X on the stack
% output = type

/Real 		2 def
/Unipotent 	1 def
/Complex 	0 def

/eigentype {
5 dict begin
/A exch def

% trace = a00 + a11
/tr A trace def

% det = a00.a11 - a01.a10
/d A det def

/Delta 
tr dup mul
4 d mul sub 
def

Delta
0 lt 
{
/type Complex def
}
if

Delta
0 gt 
{
/type Real def
}
if

Delta 
0 eq
{
A 0 get 1 get
0 eq
A 1 get 0 get
0 eq
and
{
/type Real def
} 
{
/type Unipotent def
}
ifelse

} if

type
end
} def

% - eigensolve -----------------------------------------------------

% next has input A of type Real
% leaves matrix U on stack and array of eigenvalues
% where A = UYU^{-1},Y diagonal

/real-eigensolve {
16 dict begin
/A exch def
/tr A trace def
/d A det def

/m
tr 0.5 mul 
def

/n
tr dup mul
d 4 mul
sub
sqrt
0.5 mul
def

/eigvals 
[m n add m n sub]
def

% check for diagonal matrix 
A 0 get 1 get 0 eq
A 1 get 0 get 0 eq
and
{
/x0 1 def
/y0 0 def
/x1 0 def
/y1 1 def
}
{
% set up eigenvector equations
/y0 eigvals 0 get A 0 get 0 get sub def
/x0 A 0 get 1 get def
x0 0 eq
y0 0 eq and
{
/x0 eigvals 0 get A 1 get 1 get sub def
/y0 A 1 get 0 get def
}
if
/y1 eigvals 1 get A 0 get 0 get sub def
/x1 A 0 get 1 get def
x1 0 eq
y1 0 eq and
{
/x1 eigvals 1 get A 1 get 1 get sub def
/y1 A 1 get 0 get def
}
if
} ifelse % diagonal? 

[
[x0 x1]
[y0 y1]
]

eigvals
end
} def

/unipotent-eigensolve {
16 dict begin
/A exch def
/tr A trace def
/m tr 0.5 mul def
/y0 m A 0 get 0 get sub def
/x0 A 0 get 1 get def
x0 0 eq
y0 0 eq and
{
/x0 m A 1 get 1 get sub def
/y0 A 1 get 0 get def
} if
/x0 0 eq 
{
/x1 1 def
/y1 0 def
}
{
/x1 0 def
/y1 1 def
} ifelse

[
[x0 y0]
[x1 y1]
]

m

end
} def

/complex-eigensolve {
16 dict begin
/A exch def
/d A det def
/t A trace def

/m t 0.5 mul def
/n d 4 mul t dup mul sub sqrt 0.5 mul def

[
[A 0 get 1 get 0]
[m A 0 get 0 get sub n neg]
]

[m n]
end
} def

% A 

/matrix-exponential {
8 dict begin
/A exch def
/type A 
eigentype def
type Real eq
{
(real) ==
A real-eigensolve
/eigs exch def
/X exch def
X
[
[ E eigs 0 get exp 0 ]
[ 0 E eigs 1 get exp ]
]
matrixmul
X inverse
matrixmul
} if % real
type Complex eq
{
(complex) ==
A complex-eigensolve
/eigs exch def
/X exch def
/ea E eigs 0 get exp def
/eb eigs 1 get 180 mul PI div def
X
[
[ ea eb cos mul ea eb neg sin mul ]
[ ea eb sin mul ea eb cos mul ]
]
matrixmul
X inverse matrixmul
} if % complex
type Unipotent eq {
(unipotent) ==
A unipotent-eigensolve
/eig exch def
/e E eig exp def
/X exch def
X inverse
A
matrixmul
X matrixmul
/m exch 0 get 1 get def
X 
[
[ e e m mul ]
[ 0 e]
]
matrixmul
X inverse matrixmul
} if % unipotent
end
} def

% ---------------------------------------------------------

false {

/A
[
[2 1]
[1 2]
] def

A real-eigensolve
== ==

/B
[
[2 1]
[0 2]
] def

B unipotent-eigensolve
== ==

/C 
[
[1 2]
[-2 1]
] def

C complex-eigensolve == ==

} if