Artifact dd30c71603ed7ada48d3ec454df0395edc7c0d87a75e9922c2900224639ad2d9:
- Executable file
r37/packages/normform/jordan.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 6576) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/normform/jordan.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 6576) [annotate] [blame] [check-ins using]
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;