module jordan; % Computation of the Jordan normal form of a matrix. %
load!-package 'matrix; % Otherwise setmat can fail (letmtr undefined).
load!-package 'arnum; % So we can test whether it's in use or not.
% The following functions may be useful by themselves. They can be %
% called from algebraic mode: companion, copyinto, deg_sort, %
% jordanblock, submatrix. %
% %
%**********************************************************************%
% jordan(A) computes the Jordan normal form of a square matrix A.
% First jordansymbolic is applied to A, then the symbolic zeroes of the
% characteristic polynomial are replaced by the actual zeroes. The
% zeroes of the characteristic polynomial of A are computed exactly if
% possible. The zeroes which cannot be computed exactly are
% approximated by floating point numbers.
% Specifically:
%
% - jordan(A) will return {J,P,Pinv} where J, P, and Pinv
% are such that P*J*Pinv = A.
%
% For more details of the algorithm and general workings of jordan
% see jordansymbolic.
symbolic procedure jordan(A);
begin
scalar AA,l,tmp,P,Pinv,ans,ans_mat,full_coeff_list,rule_list,
input_mode;
matrix_input_test(A);
if (car size_of_matrix(A)) neq (cadr size_of_matrix(A))
then rederr "ERROR: expecting a square matrix. ";
input_mode := get(dmode!*,'dname);
if input_mode = 'modular then rederr
"ERROR: jordan does not work with modular on. Try jordansymbolic. ";
%
% If arnum is on then we keep it on else we want
% rational on.
%
if input_mode neq 'arnum and input_mode neq 'rational
then on rational;
on combineexpt;
tmp := nest_input(A);
AA := car tmp;
full_coeff_list := cadr tmp;
l := jordansymbolicform(AA, full_coeff_list);
tmp := jordanform(l, full_coeff_list);
ans := car tmp;
P := cadr tmp;
Pinv := caddr tmp;
%
% Set up rule list for removing nests.
%
rule_list := {'co, 2,{'~, 'int}}=>'int when numberp(int);
for each elt in full_coeff_list do
<<
tmp := {'co, 2, {'~, elt}}=>elt;
rule_list := append (tmp, rule_list);
>>;
%
% Remove nests.
%
let rule_list;
ans_mat := de_nest_mat(ans);
P := de_nest_mat(P);
Pinv := de_nest_mat(Pinv);
clearrules rule_list;
%
% Return to original mode.
%
if input_mode neq 'arnum and input_mode neq 'rational then
<<
% onoff('nil,t) doesn't work so...
if input_mode = 'nil then off rational
else onoff(input_mode,t);
>>;
off combineexpt;
return {'list, ans_mat, P, Pinv};
end;
flag ('(jordan),'opfn); % So it can be used from algebraic mode.
symbolic procedure jordanform(l, full_coeff_list);
begin
scalar jj,z,zeroes,P,Pinv,x,tmp,tmp1,de_nest;
integer n,d;
P := nth(l,3);
Pinv := nth(l,4);
n := length nth(nth(l,2),1);
x := nth (nth(l,2),2);
jj := nth(l,1);
for i:=1:n do
<<
d := deg(nth(nth(nth(l,2),1),i),x);
if d>1 then
<<
%
% Determine zeroes.
%
z := nth(nth(nth(l,2),1),i);
zeroes := {};
% de-nest as solve sometimes fails with nests.
de_nest := de_nest_list(z,full_coeff_list);
tmp := algebraic solve(de_nest,x);
tmp := cdr tmp; % Remove algebraic 'list.
for j:=1:length tmp do
<<
if test_for_root_of(nth(tmp,j)) then
<<
% If poly is univariate then can solve using roots.
if length get_coeffs(de_nest) = 1 then
<<
on complex;
tmp1 := algebraic roots(z);
off complex;
>>
else
<<
on fullroots;
prin2t "***** WARNING: fullroots turned on.";
prin2t " May take a while.";
% system "sleep 3";
tmp1 := algebraic solve(de_nest,x);
off fullroots;
tmp1 := re_nest_list(tmp1,full_coeff_list);
>>;
zeroes := append(zeroes,tmp1);
zeroes := cdr zeroes;
>>
else
<<
tmp1 := algebraic solve(z,x);
tmp1 := cdr tmp1;
zeroes := append(zeroes,{nth(tmp1,j)});
>>;
>>;
%
% Substitute zeroes for symbolic names.
%
for j:=1:length zeroes do
<<
tmp := nth(zeroes,j);
tmp := caddr tmp;
jj := algebraic sub(mkid(x,off_mod_reval{'plus,
{'times,10,i},j})=tmp, jj);
>>;
for j:=1:length zeroes do
<<
tmp := nth(zeroes,j);
tmp := caddr tmp;
P := algebraic sub(mkid(x,off_mod_reval{'plus,
{'times,10,i},j})=tmp,P);
>>;
for j:=1:length zeroes do
<<
tmp := nth(zeroes,j);
tmp := caddr tmp;
Pinv := algebraic sub(mkid(x,off_mod_reval{'plus,
{'times,10,i},j})= tmp, Pinv);
>>;
>>;
>>;
return {jj,P,Pinv};
end;
symbolic procedure test_for_root_of(input);
%
% Goes through a list testing to see if there is a 'root-of
% contained within it.
%
begin
scalar tmp,copy_input,boolean,tmp1;
boolean := nil;
copy_input := input;
if atom copy_input then <<>>
else if car copy_input = 'root_of
then boolean := t
else while copy_input and boolean = nil do
<<
tmp := copy_input;
tmp := car copy_input;
if tmp = 'root_of then boolean := t
else while pairp tmp and boolean = nil do
<<
tmp1 := tmp;
if car tmp1 = 'root_of then boolean := t
else if atom tmp1 then <<>>
else while pairp tmp1 and boolean = nil do
<<
if car tmp1 = 'root_of then boolean := t
else tmp1 := car tmp1;
>>;
tmp := cdr tmp;
>>;
copy_input := cdr copy_input;
>>;
return boolean;
end;
flag ('(test_for_root_of),'boolean);
endmodule;
end;