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

module matop;

COMMENT

              #############################
              ####                     ####
              ####  MATRIX OPERATIONS  ####
              ####                     ####
              #############################


This module contains operations on dpmats, that correspond to module
operations on the corresponding images resp. cokernels.

END COMMENT;

symbolic procedure matop!=testdpmatlist l;
% Test l to be a list of dpmats embedded into a common free module.
  if null l then rederr"Empty DPMAT list"
  else begin scalar c,d;
    for each x in l do
        if not eqcar(x,'dpmat) then typerr(x,"DPMAT");
    c:=dpmat_cols car l; d:=dpmat_coldegs car l;
    for each x in cdr l do
      if not (eqn(c,dpmat_cols x) and equal(d,dpmat_coldegs x)) then
                rederr"Matrices don't match in the DPMAT list";
  end;

symbolic procedure matappend!* l;
% Appends rows of the dpmats in the dpmat list l.
  (begin scalar u,r;
      matop!=testdpmatlist l;
      cali!=degrees:=dpmat_coldegs car l;
      u:=dpmat_list car l; r:=dpmat_rows car l;
      for each y in cdr l do
        << u:=append(u, for each x in dpmat_list y collect
                        bas_newnumber(bas_nr x + r,x));
           r:=r + dpmat_rows y;
        >>;
      return dpmat_make(r,dpmat_cols car l,u,cali!=degrees,nil)
   end) where cali!=degrees:=cali!=degrees;

put('matappend,'psopfn,'matop!=matappend);
symbolic procedure matop!=matappend l;
% Append the dpmats in the list l.
  dpmat_2a matappend!* for each x in l collect dpmat_from_a reval x;

symbolic procedure mat2list!* m;
% Returns the ideal of all elements of m.
  if dpmat_cols m = 0 then m
  else (begin scalar x;
    x:=bas_renumber bas_zerodelete
        for i:=1:dpmat_rows m join
        for j:=1:dpmat_cols m collect
                bas_make(0,dpmat_element(i,j,m));
    return dpmat_make(length x,0,x,nil,
                if dpmat_cols m=1 then dpmat_gbtag m else nil)
    end) where cali!=degrees:=nil;

symbolic procedure matsum!* l;
% Returns the module sum of the dpmat list l.
  interreduce!* matappend!* l;

put('matsum,'psopfn,'matop!=matsum);
put('idealsum,'psopfn,'matop!=matsum);
symbolic procedure matop!=matsum l;
% Returns the sum of the ideals/modules in the list l.
  dpmat_2a matsum!* for each x in l collect dpmat_from_a reval x;

symbolic procedure matop!=idealprod2(a,b);
  if (dpmat_cols a > 0) or (dpmat_cols b > 0 ) then
                rederr"IDEALPROD only for ideals"
  else (begin scalar x;
    x:=bas_renumber
        for each a1 in dpmat_list a join
        for each b1 in dpmat_list b collect
            bas_make(0,dp_prod(bas_dpoly a1,bas_dpoly b1));
    return interreduce!* dpmat_make(length x,0,x,nil,nil)
    end) where cali!=degrees:=nil;

symbolic procedure idealprod!* l;
% Returns the product of the ideals in the dpmat list l.
 if null l then rederr"empty list in IDEALPROD"
 else if length l=1 then car l
 else begin scalar u;
    u:=car l;
    for each x in cdr l do u:=matop!=idealprod2(u,x);
    return u;
    end;

put('idealprod,'psopfn,'matop!=idealprod);
symbolic procedure matop!=idealprod l;
% Returns the product of the ideals in the list l.
  dpmat_2a idealprod!* for each x in l collect dpmat_from_a reval x;

symbolic procedure idealpower!*(a,n);
  if (dpmat_cols a > 0) or (not fixp n) or (n < 0) then
        rederr" Syntax : idealpower(ideal,integer)"
  else if (n=0) then dpmat_from_dpoly dp_fi 1
  else begin scalar w; w:=a;
  for i:=2:n do w:=matop!=idealprod2(w,a);
  return w;
  end;

symbolic operator idealpower;
symbolic procedure idealpower(m,l);
  if !*mode='algebraic then
        dpmat_2a idealpower!*(dpmat_from_a reval m,l)
  else idealpower!*(m,l);

symbolic procedure matop!=shiftdegs(d,n);
% Shift column degrees d n places.
   for each x in d collect ((car x + n) . cdr x);

symbolic procedure directsum!* l;
% Returns the direct sum of the modules in the dpmat list l.
  if null l then rederr"Empty DPMAT list"
  else (begin scalar r,c,u;
    for each x in l do
        if not eqcar(x,'dpmat) then typerr(x,"DPMAT")
        else if dpmat_cols x=0 then
                rederr"DIRECTSUM only for modules";
    c:=r:=0; % Actual column resp. row index.
    cali!=degrees:=nil;
    for each x in l do
       << cali!=degrees:=append(cali!=degrees,
                        matop!=shiftdegs(dpmat_coldegs x,c));
          u:=append(u, for each y in dpmat_list x collect
                bas_make(bas_nr y + r,dp_times_ei(c,bas_dpoly y)));
          r:=r + dpmat_rows x;
          c:=c + dpmat_cols x;
        >>;
    return dpmat_make(r,c,u,cali!=degrees,nil)
    end) where cali!=degrees:=cali!=degrees;

put('directsum,'psopfn,'matop!=directsum);
symbolic procedure matop!=directsum l;
% Returns the direct sum of the modules in the list l.
  dpmat_2a directsum!* for each x in l collect dpmat_from_a reval x;

symbolic operator deleteunits;
symbolic procedure deleteunits m;
  if !*noetherian then m
  else if !*mode='algebraic then dpmat_2a deleteunits!* dpmat_from_a m
  else deleteunits!* m;

symbolic procedure deleteunits!* m;
% Delete units from the base elements of the ideal m.
  if !*noetherian or (dpmat_cols m>0) then m
  else dpmat_make(dpmat_rows m,0,
        for each x in dpmat_list m collect
                bas_factorunits x,nil,dpmat_gbtag m);

symbolic procedure interreduce!* m;
  (begin scalar u;
  u:=red_interreduce dpmat_list m;
  return dpmat_make(length u, dpmat_cols m, bas_renumber u,
                cali!=degrees, dpmat_gbtag m)
  end)  where cali!=degrees:=dpmat_coldegs m;

symbolic operator interreduce;
symbolic procedure interreduce m;
% Interreduce m.
  if !*mode='algebraic then
        dpmat_2a interreduce!* dpmat_from_a reval m
  else interreduce!* m;

symbolic procedure gbasis!* m;
% Produce a minimal Groebner or standard basis of the dpmat m.
  if dpmat_gbtag m then m else car groeb_stbasis(m,t,nil,nil);

put('tangentcone,'psopfn,'matop!=tangentcone);
symbolic procedure matop!=tangentcone m;
  begin scalar c;
  intf_test m; m:=car m; intf_get m;
  if not (c:=get(m,'gbasis)) then
        put(m,'gbasis,c:=gbasis!* get(m,'basis));
  c:=tangentcone!* c;
  return dpmat_2a c;
  end;

symbolic procedure tangentcone!* m;
% Returns the tangent cone of m, provided the term order has degrees.
% m must be a gbasis.
  if null ring_degrees cali!=basering then
        rederr"tangent cone only for degree orders defined"
  else (begin scalar b;
  b:=for each x in dpmat_list m collect
    bas_make(bas_nr x,dp_tcpart bas_dpoly x);
  return dpmat_make(dpmat_rows m,
        dpmat_cols m,b,cali!=degrees,dpmat_gbtag m);
  end)  where cali!=degrees:=dpmat_coldegs m;


symbolic procedure syzygies1!* bas;
% Returns the (not yet interreduced first) syzygy module of the dpmat
% bas.
  begin
  if cali_trace() > 0 then
    << terpri(); write" Compute syzygies"; terpri() >>;
  return third groeb_stbasis(bas,nil,nil,t);
  end;

symbolic procedure syzygies!* bas;
% Returns the interreduced syzygy basis.
  interreduce!* syzygies1!* bas;

symbolic procedure normalform!*(a,b);
% Returns {a1,r,z} with a1=z*a-r*b where the rows of the dpmat a1 are
% the normalforms of the rows of the dpmat a with respect to the
% dpmat b.
   if not(eqn(dpmat_cols a,dpmat_cols b) and
        equal(dpmat_coldegs a,dpmat_coldegs b)) then
                rederr"dpmats don't match for NORMALFORM"
   else (begin scalar a1,z,u,r;
      bas_setrelations dpmat_list b;
      a1:=for each x in dpmat_list a collect
        << u:=red_redpol(dpmat_list b,x);
           z:=bas_make(bas_nr x,dp_times_ei(bas_nr x,cdr u)).z;
           car u
        >>;
      r:=bas_getrelations a1; bas_removerelations a1;
      bas_removerelations dpmat_list b; z:=reversip z;
      a1:=dpmat_make(dpmat_rows a,dpmat_cols a,a1,cali!=degrees,nil);
      cali!=degrees:=dpmat_rowdegrees b;
      r:=dpmat_make(dpmat_rows a,dpmat_rows b,bas_neworder r,
                            cali!=degrees,nil);
      cali!=degrees:=nil;
      z:=dpmat_make(dpmat_rows a,dpmat_rows a,bas_neworder z,nil,nil);
      return  {a1,r,z};
      end)  where cali!=degrees:=dpmat_coldegs b;

symbolic procedure matop_pseudomod(a,b); car mod!*(a,b);

symbolic procedure mod!*(a,b);
% Returns the normal form of the dpoly a modulo the dpmat b and the
% corresponding unit produced during pseudo division.
  (begin scalar u;
      a:=dp_neworder a; % to be on the safe side.
      u:=red_redpol(dpmat_list b,bas_make(0,a));
      return (bas_dpoly car u) . cdr u;
  end)  where cali!=degrees:=dpmat_coldegs b;

symbolic operator mod;
symbolic procedure mod(a,b);
% True normal form as s.q. also for matrices.
  if !*mode='symbolic then rederr"only for algebraic mode"
  else begin scalar u;
    b:=dpmat_from_a reval b; a:=reval a;
    if eqcar(a,'list) then
        if dpmat_cols b>0 then rederr"entries don't match for MOD"
        else a:=makelist for each x in cdr a collect
           << u:=mod!*(dp_from_a x, b);
              {'quotient,dp_2a car u,dp_2a cdr u}
           >>
    else if eqcar(a,'mat) then
        begin a:=dpmat_from_a a;
        if dpmat_cols a neq dpmat_cols b then
                rederr"entries don't match for MOD";
        a:=for each x in dpmat_list a collect mod!*(bas_dpoly x,b);
        a:='mat.
            for each x in a collect
              << u:=dp_2a cdr x;
                 for i:=1:dpmat_cols b collect
                    {'quotient,dp_2a dp_comp(i,car x),u}
              >>
        end
    else if dpmat_cols b>0 then rederr"entries don't match for MOD"
    else << u:=mod!*(dp_from_a a, b);
            a:={'quotient,dp_2a car u,dp_2a cdr u}
          >>;
    return a;
    end;

infix mod;

symbolic operator normalform;
symbolic procedure normalform(a,b);
% Compute a normal form of the rows of a with respect to b :
%   first result = third result * a + second result * b.
  if !*mode='algebraic then
  begin scalar m;
  m:= normalform!*(dpmat_from_a reval a,dpmat_from_a reval b);
  return {'list,dpmat_2a car m, dpmat_2a cadr m, dpmat_2a caddr m}
  end
  else normalform!*(a,b);

symbolic procedure eliminate!*(m,vars);
% Returns a (dpmat) basis of the elimination module of the dpmat m
% eliminating variables contained in the var. list vars.
% It sets temporary the standard elimination term order, but doesn't
% affect the ecart, and computes a Groebner basis of m.

%  if dpmat_gbtag m and eo(vars) then dpmat_sieve(m,vars,t) else

   (begin scalar c,e,bas,v;
   c:=cali!=basering; e:=ring_ecart c;
   v:=ring_names cali!=basering;
   setring!* ring_define(v,eliminationorder!*(v,vars),'revlex,e);
   cali!=degrees:=nil; % No degrees for proper result !!
   bas:=(bas_sieve(dpmat_list
            car groeb_stbasis(dpmat_neworder(m,nil),t,nil,nil), vars)
            where !*noetherian=t);
   setring!* c; cali!=degrees:=dpmat_coldegs m;
   return dpmat_make(length bas,dpmat_cols m,bas_neworder bas,
                            cali!=degrees,nil);
   end)
   where cali!=degrees:=cali!=degrees,
                cali!=basering:=cali!=basering;

symbolic operator eliminate;
symbolic procedure eliminate(m,l);
% Returns the elimination ideal/module of m with respect to the
% variables in the list l to be eliminated.
  if !*mode='algebraic then
  begin l:=reval l;
    if not eqcar(l,'list) then typerr(l,"variable list");
    m:=dpmat_from_a m; l:=cdr l;
    return dpmat_2a eliminate!*(m,l);
  end
  else eliminate!*(m,l);

symbolic procedure matintersect!* l;
  if null l then rederr"MATINTERSECT with empty list"
  else if length l=1 then car l
  else (begin scalar c,u,v,p,size;
    matop!=testdpmatlist l;
    size:=dpmat_cols car l;
    v:=for each x in l collect gensym();
    c:=cali!=basering;
    setring!* ring_sum(c,
        ring_define(v,degreeorder!* v,'lex,for each x in v collect 1));
    cali!=degrees:=mo_degneworder dpmat_coldegs car l;
    u:=for each x in pair(v,l) collect
        dpmat_times_dpoly(dp_from_a car x,dpmat_neworder(cdr x,nil));
    p:=dp_fi 1; for each x in v do p:=dp_diff(p,dp_from_a x);
    if size=0 then p:=dpmat_from_dpoly p
    else p:=dpmat_times_dpoly(p,dpmat_unit(size,cali!=degrees));
    p:=gbasis!* matsum!* (p . u);
    p:=dpmat_sieve(p,v,t);
    setring!* c;
    cali!=degrees:=dpmat_coldegs car l;
    return dpmat_neworder(p,t);
   end)
   where cali!=degrees:=cali!=degrees,
                cali!=basering:=cali!=basering;

put('matintersect,'psopfn,'matop!=matintersect);
put('idealintersect,'psopfn,'matop!=matintersect);
symbolic procedure matop!=matintersect l;
% Returns the intersection of the submodules of a fixed free module
% in the list l.
  dpmat_2a matintersect!* for each x in l collect dpmat_from_a reval x;


% ------- Submodule property and equality test --------------

put('modequalp,'psopfn,'matop!=equalp);
% Test, whether a and b are module equal.
symbolic procedure matop!=equalp u;
  if length u neq 2 then rederr"Syntax : MODEQUALP(dpmat,dpmat) "
  else begin scalar a,b;
    intf_get first u; intf_get second u;
    if null(a:=get(first u,'gbasis)) then
        put(first u,'gbasis,a:=gbasis!* get(first u,'basis));
    if null(b:=get(second u,'gbasis)) then
        put(second u,'gbasis,b:=gbasis!* get(second u,'basis));
    if modequalp!*(a,b) then return 'yes else return 'no
    end;

symbolic procedure modequalp!*(a,b);
  submodulep!*(a,b) and submodulep!*(b,a);

put('submodulep,'psopfn,'matop!=submodulep);
% Test, whether a is a submodule of b.
symbolic procedure matop!=submodulep u;
  if length u neq 2 then rederr"Syntax : SUBMODULEP(dpmat,dpmat)"
  else begin scalar a,b;
    intf_get second u;
    if null(b:=get(second u,'gbasis)) then
        put(second u,'gbasis,b:=gbasis!* get(second u,'basis));
    a:=dpmat_from_a reval first u;
    if submodulep!*(a,b) then return 'yes else return 'no
    end;

symbolic procedure submodulep!*(a,b);
  if not(dpmat_cols a=dpmat_cols b
     and equal(dpmat_coldegs a,dpmat_coldegs b)) then
    rederr"incompatible modules in SUBMODULEP"
  else (begin
    a:=for each x in dpmat_list a collect bas_dpoly x;
    return not listtest(a,b,function matop_pseudomod)
    end) where cali!=degrees:=dpmat_coldegs a;

endmodule; % matop

end;