Check-ins Using Download Hex

Artifact beff2856c4fbad705bbab87056ffe287c3558b267f970b4b3e4e4bdd0069bc25:

module jordsymb; % Computation of the Jordan normal form of a matrix.                                                              %

% The function jordansymbolic computes the Jordan normal form J of a
% matrix A, the transformation matrix P and its inverse P^(-1).
% Here symbolic names are used for the zeroes of the characteristic
% polynomial p not in K. Also a list of irreducible factors of p is
% returned.
%
% Specifically:
%
% - jordansymbolic(A) will return {J,l,P,Pinv} where J, P, and
%   Pinv are such that P*J*Pinv = A. J is calculated with symbolic 
%   names if necessary. l is {ll,x} where x is a name and ll
%   is a list of irreducible factors of p(x). If symbolic names are 
%   used then xij is a zero of ll(i).

% Global description of the algorithm:
%
% For a given n by n matrix A over a field K, we car compute the
% rational Jordan normal form R of A. Then we compute the Jordan normal
% form of R, which is also the Jordan normal form of A.
% car consider the case where R=C(p), the companion matrix of the
% monic, irreducible polynomial p=x^n+p(n-1)*x^(n-1)+..+p1*x+p0.
% If y is a zero of p then
% (y^(n-1)+p(n-1)*y^(n-2)+..+p2*y+p1,y^(n-2)+p(n-1)*y^(n-3)+..+p3*y+p2,
% ....,y^2+p(n-1)*y+p(n-2),y+p(n-1),1)
% is an eigenvector of R with eigenvalue y.
%
% Let v1     = x^(n-1)+p(n-1)*x^(n-2)+..+p2*x+p1,
%     v2     = x^(n-2)+p(n-1)*x^(n-3)+..+p3*x+p2,
%     ...
%     v(n-2) = x^2+P(n-1)*x+p(n-2),
%     v(n-1) = x+p(n-1),
%     vn     = 1.
%
% If y1,..,yn are the different zeroes of p in a splitting field of p
% over K (we asssume that p is separable, this is always true if K is a
% perfect field) we get:
%
%       inverse(V)*R*V=diag(y1,..,yn),
%
% where
%
%          [ v1(y1) v1(y2) ... v1(yn) ]
%          [ v2(y1) v2(y2) ... v2(yn) ]
%      V = [  ...    ...   ...  ...   ]
%          [  ...    ...   ...  ...   ]
%          [ vn(y1) vn(y2) ... vn(yn) ]
%
%
% One can prove that
%
%      [1 y1 ... y1^(n-1)] [v1(y1) v1(y2) ... v1(yn)]
%      [1 y2 ... y2^(n-1)] [v2(y1) v2(y2) ... v2(yn)]
%      [.................] [........................] =
%      [.................] [........................]
%      [1 yn ... yn^(n-1)] [vn(y1) vn(y2) ... vn(yn)]
%
%    = diag(diff(p,x)(y1),diff(p,x)(y2),...,diff(p,x)(yn)).
%
% If s and t are such that s*p+t*diff(p,x)=1 then we get
%
%                                           [1 y1 ... y1^(n-1) ]
%                                           [1 y2 ... y2^(n-1) ]
%    inverse(V)=diag(t(y1),t(y2),...,t(yn))*[................. ]
%                                           [................. ]
%                                           [1 yn ... yn^(n-1) ]
%
% Let Y=diag(y1,..,yn). From V^(-1)*R*V=Y it follows that
%
%                          [C(p)  I             ]
%                          [    C(p)  I         ]
%   diag(V^(-1),..,V^(-1))*[          .   .     ]*diag(V,..,V)=
%                          [            C(p)  I ]
%                          [                C(p)]
%
%          [ Y I       ]
%          [   Y I     ]
%        = [     . .   ]
%          [       Y I ]
%          [         Y ]
%
% It is now easy to see that to get our general result, we only have to
% permute diag(V,..,V) and diag(V^(-1),..,V^(-1)).

% looking_good controls formating output to print the greek character
% xi instead of lambda. At present xr does not support indexing of
% lambda but it does for all other greek letters, which is the reason
% for this switch.
%
% Only  helpful when using xr.

switch looking_good;


switch balanced_was_on;
symbolic procedure jordansymbolic(A);
  begin
    scalar AA,P,Pinv,tmp,ans_mat,ans_list,full_coeff_list,rule_list,
           output,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. ";

    if !*balanced_mod then
    <<
      off balanced_mod;
      on balanced_was_on;
    >>;

    input_mode := get(dmode!*,'dname);
    %
    % If modular or arnum are on then we keep them on else we want
    % rational on.
    %
    if input_mode neq 'modular and 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;

    tmp := jordansymbolicform(AA,full_coeff_list);
    ans_mat := car tmp;
    ans_list := cadr tmp;
    P := caddr tmp;
    Pinv := caddr rest 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_mat);
    car ans_list := append({'list},car ans_list);
    ans_list := append({'list},ans_list);
    cadr ans_list := for each elt in cadr ans_list collect reval elt;
    P := de_nest_mat(P);
    Pinv := de_nest_mat(Pinv);
    clearrules rule_list;
    %
    % Return to original mode.
    %
    if input_mode neq 'modular and 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);
       >>;
    if !*balanced_was_on then on balanced_mod;
    off combineexpt;

    output := {'list,ans_mat,ans_list,P,Pinv};
    if !*looking_good then output := looking_good(output);
    return output;
  end;

flag ('(jordansymbolic),'opfn);  %  So it can be used from 
                                 %  algebraic mode.



symbolic procedure jordansymbolicform(A,full_coeff_list);
  begin
    scalar l,R,TT,Tinv,S,Sinv,tmp,P,Pinv,invariant;

    tmp := ratjordanform(A,full_coeff_list);
    R := car tmp;
    TT := cadr tmp;
    Tinv := caddr tmp;

    tmp := ratjordan_to_jordan(R);
    l := car tmp;
    S := cadr tmp;
    Sinv := caddr tmp;

    P := off_mod_reval {'times,TT,S};
    Pinv := off_mod_reval {'times,Sinv,Tinv};

    invariant := invariant_to_jordan(nth(l,1));

    return {invariant,nth(l,2),P,Pinv};
  end;




symbolic procedure find_companion(R,rr,x);
  begin
    scalar  p;
    integer row_dim,k;

    row_dim := car size_of_matrix(R);
    k := rr+1;

    while k<=row_dim and getmat(R,k,k-1)=1 do k:=k+1;

    p := 0;
    for j:=rr:k-1 do 
    <<
      p:={'plus,p,{'times,{'minus,getmat(R,j,k-1)},{'expt,x,j-rr}}};
    >>;
    p := {'plus,p,{'expt,x,k-rr}};

    return p;
  end;




symbolic procedure find_ratjblock(R,rr,x);
  begin
    scalar p,continue;
    integer k,e,row_dim;

    row_dim := car size_of_matrix(R);
    p := reval find_companion(R,rr,x);
    e := 1;
    k:= off_mod_reval({'plus,rr,deg(p,x)});

    continue := t;
    while continue do 
    <<
      if k>row_dim then continue := nil;
      if identitymatrix(R,k-deg(p,x),k,deg(p,x)) then
      <<
        e:=e+1;
        k:=k+deg(p,x);
      >>
      else
      <<
        continue := nil;
      >>;
    >>;

    return({p,e});
  end;




symbolic procedure identitymatrix(A,i,j,m);
  %
  % Tests if the submatrix of A, starting at position (i,j) and of 
  % square size m, is an identity matrix.
  %
  begin
    integer row_dim;

    row_dim := car size_of_matrix(A);
    if i+m-1>row_dim or j+m-1>row_dim then return nil 
    else 
    <<
      if submatrix(A,{i,i+m-1},{j,j+m-1}) = make_identity(m,m) then 
         return t 
      else return nil;
    >>;
  end;

flag ('(identitymatrix),'boolean);



symbolic procedure invariant_to_jordan(invariant);
  begin
    scalar block_list;
    integer n,m;

    n := length invariant;
    block_list := {};

    for i:=1:n do
    <<
      m := length nth(nth(invariant,i),2);
      for j:=1:m do
      <<
        block_list := append(block_list,
                      {jordanblock(nth(nth(invariant,i),1),
                      nth(nth(nth(invariant,i),2),j))});
      >>;
    >>;
    return (reval {'diagi,block_list});
  end;




symbolic procedure jordanblock(const,mat_dim);
  %
  % Takes a constant (const) and an integer (mat_dim) and creates 
  % a jordan block of dimension mat_dim x mat_dim.
  %
  % Can be used independently from algebraic mode.
  %
  begin
    scalar JB;

    JB := mkmatrix(mat_dim,mat_dim);
    for i:=1:mat_dim do
    <<
      for j:=1:mat_dim do
      <<
        if i=j then 
        <<
          setmat(JB,i,j,const);
          if i<mat_dim then setmat(JB,i,j+1,1);
        >>;
      >>;
    >>;

    return JB;
  end;

flag ('(jordanblock),'opfn); %  So it can be used independently 
                             %  from algebraic mode.



switch mod_was_on;
symbolic procedure ratjordan_to_jordan(R);
  begin
    scalar prim_inv,TT,Tinv,Tinvlist,Tlist,exp_list,invariant,p,partT,
           partTinv,s1,t1,v,w,sum1,tmp,S,Sinv,x;
    integer nn,n,d;
    %
    % lambda would be better but as yet reduce can't output lambda with
    % indices (under discussion), so we use xi. If it is decided that 
    % lambda can be used with indices then change xi to lambda both 
    % here and in the functions looking_good and make_rule. 
    %
    if !*looking_good then x := 'xi
    else x := 'lambda;

    prim_inv := ratjordan_to_priminv(R,x);

    invariant := {};
    Tlist := {};
    Tinvlist := {};

    nn := length prim_inv; 
    for i:=1:nn do
    <<
      p := nth(nth(prim_inv,i),1);
      exp_list := nth(nth(prim_inv,i),2);
      d := off_mod_reval(deg(p,x));
      if d=1 then
      <<
        invariant := append(invariant,{{reval {'minus,coeffn(p,x,0)},
                            exp_list}});
      >>
      else
      <<
        for j:=1:d do
        <<
          invariant := append(invariant,{{mkid(x,off_mod_reval{'plus,
                              {'times,10,i},j}),exp_list}});
        >>;
      >>;
      %
      % Compute eigenvector of C(p) with eigenvalue x.
      %
      v := mkvect(d);
      putv(v,d,1);
      for j:=d-1 step -1 until 1 do
      <<
        tmp := 0;
        sum1 := 0;
        for k:=j:d-1 do
        <<
          tmp := reval{'times,coeffn(p,x,k),{'expt,x,k-j}};
          sum1 := reval{'plus,sum1,tmp};
        >>;
        putv(v,j,reval {'plus,sum1,{'expt,x,d-j}});
      >>;

      sum1 := 0;
      for j:=1:length exp_list do
      <<
        tmp := reval nth(exp_list,j);
        sum1 := reval {'plus,sum1,tmp};
      >>;
      n := sum1;

      partT := mkmatrix(n*d,n);
      for j:=1:n do
      <<
        for k:=1:d do
        <<
          setmat(partT,(j-1)*d+k,j,getv(v,k));
        >>;
      >>;

      TT := mkmatrix(n*d,n*d);
      if d=1 then
      <<
        %
        % off modular else the mkid number part will be calculated
        % in modular mode.
        %
        if !*modular then << off modular; on mod_was_on; >>;
        TT := copyinto(algebraic (sub({x=-coeffn(p,x,0)},partT)),
                       TT,1,1);
        if !*mod_was_on then << on modular; off mod_was_on; >>;
      >>
      else for j:=1:d do
      <<
        %
        % off modular else the mkid number part will be calculated 
        % in modular mode.
        %
        if !*modular then << off modular; on mod_was_on; >>;
        TT := copyinto(algebraic sub(x=mkid(x,off_mod_reval{'plus,
                       {'times,10,i},j}),partT),TT,1,(j-1)*n+1);
        if !*mod_was_on then << on modular; off mod_was_on; >>;
      >>;
      Tlist := append(Tlist,{TT});

      tmp := algebraic df(p,x);
      tmp := calc_exgcd(p,tmp,x);
      s1 := cadr tmp;
      t1 := caddr tmp;
      w := mkvect(d);
      putv(w,1,t1);

      for j:=2:d do
      <<
        putv(w,j,get_rem(reval{'times,x,getv(w,j-1)},p));
      >>;
  
      partTinv := mkmatrix(n,n*d);
      for j:=1:n do
      <<
        for k:=1:d do
        <<
          setmat(partTinv,j,(j-1)*d+k,getv(w,k));
        >>;
      >>;

      Tinv := mkmatrix(n*d,n*d);
      if d=1 then
      <<
        %
        % off modular else the mkid number part will be calculated
        % in modular mode.
        %
        if !*modular then << off modular; on mod_was_on; >>;
        Tinv := reval copyinto(algebraic sub(x=-coeffn(p,x,0),partTinv)
                               ,Tinv,1,1);
        if !*mod_was_on then << on modular; off mod_was_on; >>;
      >>
      else for j:=1:d do
      <<
        %
        % off modular else the mkid number part will be calculated
        % in modular mode.
        %
        if !*modular then << off modular;  on mod_was_on; >>;
        Tinv := reval copyinto(algebraic sub(x=mkid(x,off_mod_reval
                {'plus,{'times,10,i},j}),partTinv),Tinv,(j-1)*n+1,1);
        if !*mod_was_on then << on modular; off mod_was_on; >>;
      >>;
      Tinvlist := append(Tinvlist,{Tinv});

    >>;

    S := reval{'diagi,Tlist};
    Sinv := reval{'diagi,Tinvlist};
    tmp := {for i:=1:nn collect nth(nth(prim_inv ,i),1)};
    tmp := append(tmp,{x});
    tmp := append({invariant},{tmp});

    return {tmp,S,Sinv};
  end;




symbolic procedure ratjordan_to_priminv(R,x);
  %
  % ratjordan_to_priminv(R,x) computes the primary invariant of a matrix
  % R which is in rational Jordan normal form.
  %
  begin
    scalar p,plist,exp_list,l,prim_inv;
    integer n,rr,ii,nn;

    n := car size_of_matrix(R);
    rr := 1;

    plist := {};
    while rr<=n do
    <<
      l := find_ratjblock(R,rr,x);
      plist := append(plist,{l});
      rr := off_mod_reval({'plus,rr,{'times,nth(l,2),deg(nth(l,1),x)}});
    >>;

    p := reval nth(nth(plist,1),1);
    exp_list := {nth(nth(plist,1),2)};
    prim_inv := {};
    nn := length plist;
    ii := 2;
    while ii<=nn do
    <<
      if reval nth(nth(plist,ii),1) = p then
      <<
        exp_list := append(exp_list,{nth(nth(plist,ii),2)})
      >>
      else
      <<
        prim_inv := append(prim_inv,{{p,exp_list}});
        p := reval nth(nth(plist,ii),1);
        exp_list := {nth(nth(plist,ii),2)};
      >>;
      ii := ii+1;
    >>;
    prim_inv := append(prim_inv,{{p,exp_list}});

    return prim_inv;
  end;




symbolic procedure submatrix(A,row_list,col_list);
  %
  % Creates the submatrix of A from rows  p to q (row_list = {p,q})
  % and columns r to s (col_list = {r,s}).
  %
  % Can be used independently from algebraic mode.
  %
  begin
    scalar AA;
    integer row_dim,col_dim,car_row,last_row,car_col,last_col,
            A_row_dim,A_col_dim;

    matrix_input_test(A);

    %  If algebraic input remove 'list.
    if car row_list = 'list then row_list := cdr row_list;
    if car col_list = 'list then col_list := cdr col_list;

    car_row := car row_list;
    last_row := cadr row_list;
    row_dim := last_row - car_row + 1;
    car_col := car col_list;
    last_col := cadr col_list;
    col_dim := last_col - car_col + 1;

    A_row_dim := car size_of_matrix(A);
    A_col_dim := cadr size_of_matrix(A);
    if car_row = 0 or last_row = 0 then rederr
     {"0 is out of range for ",A,". The car row is labelled 1."};
    if car_col = 0 or last_col = 0 then rederr
     {"0 is out of range for",A,". The car column is labelled 1."};
    if car_row > A_row_dim then rederr 
     {A,"doesn't have",car_row,"rows."};
    if last_row > A_row_dim then rederr 
     {A,"doesn't have",last_row,"rows."};
    if car_col > A_col_dim then rederr 
     {A,"doesn't have",car_col,"columns."};
    if last_col > A_col_dim then rederr 
     {A,"doesn't have",last_col,"columns."};

    AA := mkmatrix(row_dim,col_dim);
    for i:=1:row_dim do
    <<
      for j:=1:col_dim do
      <<
        setmat(AA,i,j,getmat(A,i+car_row-1,j+car_col-1));
      >>;
    >>;

    return AA;
  end;

flag ('(submatrix),'opfn);  %  So it can be used independently 
                            %  from algebraic mode.



symbolic procedure looking_good(output);
  %
  % Converts output for correct display of indices with greek 
  % font. Only used when switch looking_good is on, which is only on 
  % when using xr.
  %
  % xiab => xi(a,b) is used instead of xiab => xi(ab) to reduce problems
  % when using modular arithmetic. In mod 17 (for example) xi(21) will 
  % get converted to xi(2,1) but the alternative would give xi(4) - 
  % wrong! Unfortunately the alternative probably looks a bit nicer but 
  % there you go. If the modulus is <= 7 then this may give problems, 
  % eg: xi55 in mod 5 will give xi(0,0). These cases will be very rare 
  % but if they occur turn OFF looking_good.
  % 
  begin
    scalar rule_list;

    algebraic operator xi;
    algebraic print_indexed(xi);
    %
    % Create rule list to convert xiab => xi(a,b) for a,b:=1:9.
    % 
    rule_list := make_rule();

    let rule_list;
    output := reval output;
    clearrules rule_list;

    return output;
  end;




symbolic procedure make_rule();
  begin
    scalar rule_list,tmp,tmp1;

    rule_list := {};
    for i:=1:9 do
      <<
        for j:=1:9 do
        <<
          tmp1 := reval mkid('xi,i);
          tmp1 := reval mkid(tmp1,j);
          tmp := {'replaceby,tmp1,{'xi,i,j}};
          rule_list := append(rule_list,{tmp});
        >>;
      >>;
      rule_list := append({'list},rule_list);
      return rule_list;
  end;

endmodule;  

end;
REDUCE Historical
[ REDUCE Sourceforge Project | Historical SVN Repository | GitHub Mirror | SourceHut Mirror | NotABug Mirror | Chisel Mirror | Chisel RSS ]