reduce-historical/r37/packages/cali/dpmat.red

module dpmat;

COMMENT

                 #####################
                 ###               ###
                 ###    MATRICES   ###
                 ###               ###
                 #####################

This module introduces special dpoly matrices with its own matrix
syntax.

Informal syntax :

 <matrix> ::= list('DPMAT,#rows,#cols,baslist,column_degrees,gb-tag)

Dpmat's are the central data structure exploited in the modules of
this package. Each such matrix describes a map

                   f : R^rows --> R^cols,

that gives rise for the definition of two modules,

                im f = the submodule of R^cols generated by the rows
                       of the matrix

        and  coker f = R^cols/im f.

Conceptually dpmat's are identified with im f.

END COMMENT;

% ------------- Reference operators ----------------

symbolic procedure dpmat_rows m; cadr m;
symbolic procedure dpmat_cols m; caddr m;
symbolic procedure dpmat_list m; cadddr m;
symbolic procedure dpmat_coldegs m; nth(m,5);
symbolic procedure dpmat_gbtag m; nth(m,6);

% ------------- Elementary operations --------------

symbolic procedure dpmat_rowdegrees m;
% Returns the row degrees of the dpmat m as an assoc. list.
  (for each x in dpmat_list m join
        if (bas_nr x > 0) and bas_dpoly x then
                {(bas_nr x).(mo_getdegree(dp_lmon bas_dpoly x,l))})
  where l=dpmat_coldegs m;

symbolic procedure dpmat_make(r,c,bas,degs,gbtag);
  list('dpmat,r,c,bas,degs,gbtag);

symbolic procedure dpmat_element(r,c,mmat);
% Returns mmat[r,c].
   dp_neworder
        dp_comp(c, bas_dpoly bas_getelement(r,dpmat_list mmat));

symbolic procedure dpmat_print m; mathprint dpmat_2a m;

symbolic procedure getleadterms!* m;
% Returns the dpmat with the leading terms of m.
  (begin scalar b;
  b:=for each x in dpmat_list m collect
    bas_make(bas_nr x,list(car bas_dpoly x));
  return dpmat_make(dpmat_rows m,dpmat_cols m,b,cali!=degrees,t);
  end) where cali!=degrees:=dpmat_coldegs m;

% -------- Symbolic mode file transfer --------------

symbolic procedure savemat!*(m,name);
% Save the dpmat m under the name <name>.
  begin scalar nat,c;
  if not (stringp name or idp name) then typerr(name,"file name");
  if not eqcar(m,'dpmat) then typerr(m,"dpmat");
  nat:=!*nat; !*nat:=nil;
  write"Saving as ",name;
  out name$
  write"algebraic(setring "$

  % mathprint prints lists without terminator, but matrices with
  % terminator.

  mathprint ring_2a cali!=basering$ write")$"$
  write"algebraic(<<basis :="$ dpmat_print m$
  if dpmat_cols m=0 then write"$"$ write">>)$"$
  if (c:=dpmat_coldegs m) then
    << write"algebraic(degrees:="$
    mathprint moid_2a for each x in c collect cdr x$ write")$"$
    >>;
  write"end$"$ terpri()$
  shut name; terpri(); !*nat:=nat;
  end;

symbolic procedure initmat!* name;
% Initialize a dpmat from <name>.
  if not (stringp name or idp name) then typerr(name,"file name")
  else begin scalar m,c; integer i;
     write"Initializing ",name; terpri();
     in name$ m:=reval 'basis; cali!=degrees:=nil;
     if eqcar(m,'list) then
        << m:=bas_from_a m; m:=dpmat_make(length m,0,m,nil,nil)>>
     else if eqcar(m,'mat) then
        << c:=moid_from_a reval 'degrees; i:=0;
           cali!=degrees:=for each x in c collect <<i:=i+1; i . x >>;
           m:=dpmat_from_a m;
        >>
     else typerr(m,"basis or matrix");
     dpmat_print m;
     return m;
     end;

% ---- Algebraic mode file transfer ---------

symbolic operator savemat;
symbolic procedure savemat(m,name);
  if !*mode='algebraic then savemat!*(dpmat_from_a m,name)
  else savemat!*(m,name);

symbolic operator initmat;
symbolic procedure initmat name;
  if !*mode='algebraic then dpmat_2a initmat!* name
  else initmat!* name;

% --------------- Arithmetics for dpmat's ----------------------

symbolic procedure dpmat!=dpsubst(a,b);
% Substitute in the dpoly a each e_i by b_i from the base list b.
   begin scalar v;
   for each x in b do
        v:=dp_sum(v,dp_prod(dp_comp(bas_nr x,a),bas_dpoly x));
   return v;
   end;

symbolic procedure dpmat_mult(a,b);
% Returns a * b.
  if not eqn(dpmat_cols a,dpmat_rows b) then
        rerror('dpmat,1," matrices don't match for MATMULT")
  else dpmat_make( dpmat_rows a, dpmat_cols b,
        for each x in dpmat_list a collect
                bas_make(bas_nr x,
                        dpmat!=dpsubst(bas_dpoly x,dpmat_list b)),
        cali!=degrees,nil)
        where cali!=degrees:=dpmat_coldegs b;

symbolic procedure dpmat_times_dpoly(f,m);
% Returns f * m for the dpoly f and the dpmat m.
   dpmat_make(dpmat_rows m,dpmat_cols m,
        for each x in dpmat_list m collect
            bas_make1(bas_nr x, dp_prod(f,bas_dpoly x),
                        dp_prod(f,bas_rep x)),
        cali!=degrees,nil) where cali!=degrees:=dpmat_coldegs m;

symbolic procedure dpmat_neg a;
% Returns - a.
   dpmat_make(
        dpmat_rows a,
        dpmat_cols a,
        for each x in dpmat_list a collect
            bas_make1(bas_nr x,dp_neg bas_dpoly x, dp_neg bas_rep x),
        cali!=degrees,dpmat_gbtag a)
        where cali!=degrees:=dpmat_coldegs a;

symbolic procedure dpmat_diff(a,b);
% Returns a - b.
  dpmat_sum(a,dpmat_neg b);

symbolic procedure dpmat_sum(a,b);
% Returns a + b.
  if not (eqn(dpmat_rows a,dpmat_rows b)
                and eqn(dpmat_cols a, dpmat_cols b)
                and equal(dpmat_coldegs a,dpmat_coldegs b)) then
        rerror('dpmat,2,"matrices don't match for MATSUM")
  else (begin scalar u,v,w;
        u:=dpmat_list a; v:=dpmat_list b;
        w:=for i:=1:dpmat_rows a collect
            (bas_make1(i,dp_sum(bas_dpoly x,bas_dpoly y),
                            dp_sum(bas_rep x,bas_rep y))
                        where y= bas_getelement(i,v),
                              x= bas_getelement(i,u));
        return dpmat_make(dpmat_rows a,dpmat_cols a,w,cali!=degrees,
                          nil);
        end) where cali!=degrees:=dpmat_coldegs a;

symbolic procedure dpmat_from_dpoly p;
  if null p then dpmat_make(0,0,nil,nil,t)
  else dpmat_make(1,0,list bas_make(1,p),nil,t);

symbolic procedure dpmat_unit(n,degs);
% Returns the unit dpmat of size n.
  dpmat_make(n,n, for i:=1:n collect bas_make(i,dp_from_ei i),degs,t);

symbolic procedure dpmat_unitideal!? m;
  (dpmat_cols m = 0) and null matop_pseudomod(dp_fi 1,m);

symbolic procedure dpmat_transpose m;
% Returns transposed m with consistent column degrees.
  if (dpmat_cols m = 0) then dpmat!=transpose ideal2mat!* m
  else dpmat!=transpose m;

symbolic procedure dpmat!=transpose m;
  (begin scalar b,p,q;
     cali!=degrees:=
        for each x in dpmat_rowdegrees m collect
                (car x).(mo_neg cdr x);
     for i:=1:dpmat_cols m do
       << p:=nil;
          for j:=1:dpmat_rows m do
            << q:=dpmat_element(j,i,m);
               if q then p:=dp_sum(p,dp_times_ei(j,q))
            >>;
          if p then b:=bas_make(i,p) . b;
       >>;
   return dpmat_make(dpmat_cols m,dpmat_rows m,reverse b,
                        cali!=degrees,nil);
   end) where cali!=degrees:=cali!=degrees;

symbolic procedure ideal2mat!* u;
% Returns u as column vector if dpmat_cols u = 0.
  if dpmat_cols u neq 0 then
            rerror('dpmat,4,"IDEAL2MAT only for ideal bases")
  else dpmat_make(dpmat_rows u,1,
                 for each x in dpmat_list u collect
                   bas_make(bas_nr x,dp_times_ei(1,bas_dpoly x)),
                nil,dpmat_gbtag u) where cali!=degrees:=nil;

symbolic procedure dpmat_renumber old;
% Renumber dpmat_list old.
% Returns (new . change) with new = change * old.
  if null dpmat_list old then (old . dpmat_unit(dpmat_rows old,nil))
  else (begin scalar i,u,v,w;
      cali!=degrees:=dpmat_rowdegrees old;
      i:=0; u:=dpmat_list old;
      while u do
        <<i:=i+1; v:=bas_newnumber(i,car u) . v;
        w:=bas_make(i,dp_from_ei bas_nr car u) . w  ; u:=cdr u>>;
      return dpmat_make(i,dpmat_cols old,
                        reverse v,dpmat_coldegs old,dpmat_gbtag old) .
             dpmat_make(i,dpmat_rows old,reverse w,cali!=degrees,t);
      end) where cali!=degrees:=cali!=degrees;

symbolic procedure mathomogenize!*(m,var);
% Returns m with homogenized rows using the var. name var.
  dpmat_make(dpmat_rows m, dpmat_cols m,
        bas_homogenize(dpmat_list m,var), cali!=degrees,nil)
  where cali!=degrees:=dpmat_coldegs m;

symbolic operator mathomogenize;
symbolic procedure mathomogenize(m,v);
% Returns the homogenized matrix of m with respect to the variable v.
  if !*mode='algebraic then
        dpmat_2a mathomogenize!*(dpmat_from_a reval m,v)
  else matdehomogenize!*(m,v);

symbolic procedure matdehomogenize!*(m,var);
% Returns m with var. name var set equal to one.
  dpmat_make(dpmat_rows m, dpmat_cols m,
        bas_dehomogenize(dpmat_list m,var), cali!=degrees,nil)
  where cali!=degrees:=dpmat_coldegs m;

symbolic procedure dpmat_sieve(m,vars,gbtag);
% Apply bas_sieve to dpmat_list m. The gbtag slot allows to set the
% gbtag of the result.
  dpmat_make(length x,dpmat_cols m,x,cali!=degrees,gbtag)
        where x=bas_sieve(dpmat_list m,vars)
  where cali!=degrees:=dpmat_coldegs m;

symbolic procedure dpmat_neworder(m,gbtag);
% Apply bas_neworder to dpmat_list m with current cali!=degrees.
% The gbtag sets the gbtag part of the result.
   dpmat_make(dpmat_rows m,dpmat_cols m,
        bas_neworder dpmat_list m,cali!=degrees,gbtag);

symbolic procedure dpmat_zero!? m;
% Test whether m is a zero dpmat.
   bas_zero!? dpmat_list m;

symbolic procedure dpmat_project(m,k);
% Project the dpmat m onto its first k components.
  dpmat_make(dpmat_rows m,k,
        for each x in dpmat_list m collect
                bas_make(bas_nr x,dp_project(bas_dpoly x,k)),
        dpmat_coldegs m,nil);

% ---------- Interface to algebraic mode

symbolic procedure dpmat_2a m;
% Convert the dpmat m to a matrix (c>0) or a polynomial list (c=0) in
% algebraic (pseudo)prefix form.
  if dpmat_cols m=0 then bas_2a dpmat_list m
  else   'mat .
        if dpmat_rows m=0 then list for j:=1:dpmat_cols m collect 0
        else  for i:=1:dpmat_rows m collect
                for j:=1:dpmat_cols m collect
                    dp_2a dpmat_element(i,j,m);

symbolic procedure dpmat_from_a m;
% Convert an algebraic polynomial list or matrix expression into a
% dpmat with respect to the current setting of cali!=degrees.
  if eqcar(m,'mat) then
    begin integer i; scalar u,p; m:=cdr m;
    for each x in m do
      << i:=1; p:=nil;
         for each y in x do
           << p:=dp_sum(p,dp_times_ei(i,dp_from_a reval y)); i:=i+1 >>;
         u:=bas_make(0,p).u
      >>;
    return dpmat_make(length m,length car m,
                bas_renumber reversip u, cali!=degrees,nil);
    end
  else if eqcar(m,'list) then
    ((begin scalar x;  x:=bas_from_a reval m;
    return dpmat_make(length x,0,x,nil,nil)
    end) where cali!=degrees:=nil)
  else typerr(m,"polynomial list or matrix");

% ---- Substitution in dpmats --------------

symbolic procedure dpmat_sub(a,m);
% a=list of (var . alg. prefix form) to be substituted into the dpmat
% m.
   dpmat_from_a subeval1(a,dpmat_2a m)
  where cali!=degrees:=dpmat_coldegs m;

% ------------- Determinant ------------------------

symbolic procedure dpmat_det m;
% Returns the determinant of the dpmat m.
  if dpmat_rows m neq dpmat_cols m then rederr "non-square matrix"
  else dp_from_a prepf numr detq matsm dpmat_2a m;

endmodule; % dpmat

end;