Check-ins Using Download Hex

Artifact 1258ee1d24e838c49d33db174c2a9f644f17163d8d4adfc0a98fa68445f4be7f:

module xideal;

%
%                               XIDEAL
%
% Tools for calculations with ideals of polynomials in exterior algebra.
% Uses Groebner basis algorithms described in D Hartley and P A Tuckey,
% "A direct characterisation of Groebner bases in Clifford and Grassmann
% algebras", Preprint MPI-Ph/93-96 1993.
%
% Authors:      David Hartley           Philip A Tuckey
%               GMD, Institute I1       Max Planck Institute for Physics
%               Schloss Birlinghoven    Foehringer Ring 6
%               D-53757 St Augustin     D-80805 Munich
%               Germany                 Germany
%
%               email:  David.Hartley@gmd.de
%                       pht@iws170.mppmu.mpg.de
%
% Created:      1/5/92          as idl.red
%
% Modified:     5/8/92          Renamed ideal.red
%                               Chose xideal30 as algorithm
%                               Renamed mod -> xmodulo, ideal -> xideal
%                               Added xmodulop
%                               Added subform
%                               Added autoreduction (xfullreduce)
%               4/3/94          Renamed xideal.red
%                               Compiles independently
%                               Converted right reduction and spolys to
%                                  left
%                               Added graded lexicographical ordering
%                               Enabled non-graded ideals
%                               Fixed trivial ideal bug
%                               Removed subform
%                               Renamed xtrace -> xstats
%
%
% Algebraic mode entry points
%
% xideal(S) or xideal(S,r)
%  calculates an exterior Groebner basis for the list of generator S,
%  optionally up to degree r.
%
% xmodulo(F,S) or F xmodulo S
% reduces F with respect to an exterior Groebner basis for the list of
% generators S.  In the first format, F may be either a single exterior
% form, or a list of forms.  If this routine is used more than once,
% and S does not change between calls, then the Groebner basis is not
% recalculated.

% xmodulop(F,S) or F xmodulop S
% reduces F with respect to the set of exterior polynomials S, which is
% not necessarily a Groebner basis.  In the first format, F may be
% either a single exterior form, or a list of forms.  This routine can
% be used in conjunction with xideal to produce the same effect as
% xmodulo:  F xmodulo S = F xmodulop xideal(S,degree F).

% Switches
%
% xstats        - Produces counting and timing information (default OFF)
% xfullreduce   - Allows reduced Groebner bases to be calculated
%                    (default ON)

% ======================================================================


% ----------------------------------------------------------------------

create!-package('(xideal xutils),'(contrib xideal));

fluid '(!*xstats !*xfullreduce !*xdebug
        xtruncate!* xideal!* xbasis!* xdegree!*
        xidealtime!* xspolytime!* xwedgestime!*
        xreducetime!* xautoreducetime!*
        xwedges!* newxwedges!* xspolys!* newxspolys!*);

!*xdebug := nil;        % non-nil prints algorithm trace
xtruncate!* := nil;     % non-nil gives truncation degree for GB

xbasis!* := nil;        % last input set to xmodulo
xideal!* := nil;        % last output from xmodulo
xdegree!* := nil;       % largest recent degree to xmodulo

switch 'xstats;
!*xstats := nil;

switch 'xfullreduce;
!*xfullreduce := t;

infix xmodulo;
precedence xmodulo,freeof;

infix xmodulop;
precedence xmodulop,freeof;


% ----------------------------------------------------------------------

% These definitions copied from E. Schruefer's EXCALC for compilation.

global '(dimex!*);

smacro procedure ldpf u;
   %selector for leading standard form in patitioned sf;
   caar u;

smacro procedure !*k2pf u;
   u .* (1 ./ 1) .+ nil;

smacro procedure negpf u;
   multpfsq(u,(-1) ./ 1);

% ----------------------------------------------------------------------


symbolic operator xmodulo;
symbolic procedure xmodulo(x,l);
% x is prefix form or algebraic list, l is algebraic list, result is
% (possibly algebraic list of) !*sq prefix
% if generators haven't changed, use last ideal (speeds up loops)
begin scalar z;
 l := !*al2sl l;
 if atom x or car x neq 'list then
   z := x
 else   % extract maximum degree from list
   <<z := 0;
     x := 'list . !*al2sl x;
     for each y in cdr x do if xdegree y > xdegree z then z := y>>;
 xtruncate!* := xdegree!*; % just in case it was used in-between
 if (l neq xbasis!*) or xdegreecheck z then
   <<xideal!* := xideal0(l,xdegree z);
     xdegree!* := xdegree z;
     xbasis!* := l>>;
 z := if atom x or car x neq 'list then list x else cdr x;
 z := for each y in z collect
        mk!*sq !*pf2sq xreduce(xreorder partitop y,xideal!*);
 return if atom x or car x neq 'list then car z
         else purge strip0('list . z);
end;


% ----------------------------------------------------------------------


symbolic operator xmodulop;
symbolic procedure xmodulop(x,l);
% x is prefix form or algebraic list, l is algebraic list, result is
% (possibly algebraic list of) !*sq prefix
begin
 l := for each q in !*al2sl l collect xreorder partitop q;
 if atom x or car x neq 'list then
   return mk!*sq !*pf2sq xreduce(xreorder partitop x,l);
 return purge strip0('list . for each y in !*al2sl x collect
                         mk!*sq !*pf2sq xreduce(xreorder partitop y,l));
end;


% ----------------------------------------------------------------------


put('xideal,'psopfn,'xidealeval);

symbolic procedure xidealeval u;
% u is a list with one or two items: an algebraic list of p-forms
% and an integer
begin scalar l,r,n;
 n := length u;
 if n < 1 or n > 2 then rederr "Wrong number of arguments to XIDEAL";
 l := !*al2sl reval car u;
 if n = 2 then r := reval cadr u;
return 'list . for each f in xideal0(l,r) collect mk!*sq !*pf2sq f;
end;


% ----------------------------------------------------------------------


symbolic procedure xideal0(p,r);
% p is list of prefix, r is an integer or nil, result is list of pf's

xidealpf(for each f in p collect partitop f,r);


% ----------------------------------------------------------------------


symbolic procedure xidealpf(p,r);
% p is list of pf's, r is an integer or nil, result is list of pf's
begin scalar p1,p2,d,nongraded,vars;
 rmsubs();
 if !*xstats then xstatsinitialise();
 xtruncate!* := r;
 p := nil . p;
 while (p := cdr p) do
   begin
    vars := xvarcheck(car p,vars);
    d := xhomogeneous car p;
    if null d then <<nongraded := t; xtruncate!* := nil>>;
    if d = 0 then  % quick escape, trivial ideal
      <<p := p1 := list partitop 1; p2 := nil>>
    else if d = 1 then p1 := xreorder car p . p1
    else p2 := xreorder car p . p2;
  end;
 p1 := xautoreduce p1;
 for each s in p2 do
   if (s := xreduce(s,append(p1,p))) and not xdegreecheck ldpf s then
     begin
      p := xidealpf0(p,s);
      if !*xstats
        then <<prin2 "GB size  : "; prin2t(length p + length p1)>>;
      if !*xfullreduce then p := xautoreduce p;
      if !*xstats and !*xfullreduce then
        <<prin2 "RGB size : "; prin2t(length p + length p1)>>;
     end;
 p := append(p1,p);
 if nongraded and !*xfullreduce and p1 then p := xautoreduce p;
 if !*xstats then xstatsprint();
 return p;
end;


% ----------------------------------------------------------------------


symbolic procedure xstatsinitialise;
begin
 xidealtime!* := time(); xspolytime!* := 0; xwedgestime!* := 0;
 xreducetime!* := 0; xautoreducetime!* := 0;
 xwedges!* := 0; newxwedges!* := 0;
 xspolys!* := 0; newxspolys!* := 0;
end;


% ----------------------------------------------------------------------


symbolic procedure xstatsprint;
begin
 xidealtime!* := time() - xidealtime!*;
 prin2 "Wedges generated     : "; prin2t xwedges!*;
 prin2 "Wedges surviving     : "; prin2t newxwedges!*;
 prin2 "Survival rate        : ";
 if xwedges!* = 0 then prin2 0 else prin2((100*newxwedges!*)/xwedges!*);
 prin2t "%";
 prin2 "Spolys generated     : "; prin2t xspolys!*;
 prin2 "Spolys surviving     : "; prin2t newxspolys!*;
 prin2 "Survival rate        : ";
 if xspolys!* = 0 then prin2 0 else prin2((100*newxspolys!*)/xspolys!*);
 prin2t "%";
 prin2 "Total computing time : "; prin2 xidealtime!*; prin2t " ms";
 prin2 "Spoly time           : "; prin2 xspolytime!*; prin2t " ms";
 prin2 "Wedges time          : "; prin2 xwedgestime!*; prin2t " ms";
 prin2 "Reduction time       : "; prin2 xreducetime!*; prin2t " ms";
 prin2 "Auto-reduction time  : "; prin2 xautoreducetime!*; prin2t " ms";
end;


% ----------------------------------------------------------------------


symbolic procedure xidealpf0(p,s);
% p is list of pf's, s is pf, result is list of pf's
% p is a GB, result is a GB containing p and s.
begin scalar h,b,y;
 if !*xstats then
   <<newxspolys!* := newxspolys!* - 1;
     princ '#; if !*xdebug then prin2t s>>;
 while s do
   begin
    if !*xstats then newxspolys!* := newxspolys!* + 1;
    y := xwedges(s,append(p,h));
    if ldpf car y = 1 then % quick escape, trivial ideal
      return <<s := nil; h := nil; p := list partitop 1;>>;
    for each f in y do
      <<b := append(b,
             append(for each j in p collect list(j,f),
                    for each j in h collect list(j,f)));
        h := append(h,list f)>>;
    s := nil;
    while b and null s do
      <<s := xreduce(xspoly(car b),append(p,h));
        if !*xstats and s then
          <<princ "$";
            if !*xdebug then
              <<prin2t caar b; prin2 "*"; prin2t cadar b;
                prin2 " --> "; prin2t s>> >>;
        b := cdr b>>;
   end;
 if !*xstats then prin2t "";
 return append(p,h);
end;


% ----------------------------------------------------------------------


symbolic procedure xspoly p;
% p is a list of two reordered pf's, result is reordered pf
begin scalar u,v,lu,lv,cu,cv,z,t0;
 if !*xstats then t0 := time();
 u := car p; v := cadr p;
 if null u or null v then return nil;
 lu := ldpffax u;
 lv := ldpffax v;
 cu := setdiff(lu,lv);
 cv := setdiff(lv,lu);
 if (cu = lu) % cv = lv as well in this case
    or xdegreecheck('wedge . append(lu,cv)) then return nil;
 if null cu then cu := list 1;  if null cv then cv := list 1;
 u := partitwedge list('wedge . cv,mk!*sq !*pf2sq u);
 v := partitwedge list('wedge . cu,mk!*sq !*pf2sq v);
 z := negsq quotsq(lc u,lc v);
 z := xreorder addpf(u,multpfs(1 .* z .+ nil,v));
 if !*xstats then xspolys!* := xspolys!* + 1;
 if !*xstats then xspolytime!* := xspolytime!* + time() - t0;
 return z;
end;


% ----------------------------------------------------------------------


symbolic procedure xwedges(s,p);
% s is reordered pf, p is a list of reordered pf's,
% result is a list of reordered pf's (first element is always s, unless
% trivial ideal detected)
begin scalar t0,y,h,q,l;
 if xdegreecheck ldpf s then
   rederr "Major problem -- please contact XIDEAL authors";
 if d and d leq 1 where d = xhomogeneous s then return list s;
 if !*xstats then t0 := time();
 if xtruncate!* then xtruncate!* := xtruncate!* - 1;
 y := list s;
 h := p;
 while y do
   begin
    s := car y;
    y := cdr y;
    h := append(h,list s);
    if xdegreecheck ldpf s or null red s then return;
    l := nil . ldpffax s;
    while (l := cdr l) do
      begin
       if !*xstats then xwedges!* := xwedges!* + 1;
       q := xreorder wedgepf(!*k2pf car l,s);
       if null (q := xreduce(q,append(h,y))) then return;
       if !*xstats then
         <<newxwedges!* := newxwedges!* + 1;
           princ "^";
           if !*xdebug then
             <<prin2 car l; prin2 " --> "; prin2t q>> >>;
       if ldpf q = 1 then % quick escape, trivial ideal
         return <<l := list nil; y := nil; % terminate loops
                  h := append(p,list partitop 1)>>;
       y := append(y, list q);
      end;
   end;
 if xtruncate!* then xtruncate!* := xtruncate!* + 1;
 if !*xstats then xwedgestime!* := xwedgestime!* + time() - t0;
 return setdiff(h,p);
end;


% ----------------------------------------------------------------------


symbolic procedure xautoreduce0 p;
% p is a list of pf's, result is a list of pf's
% symbolic mode entry point to xautoreduce for unordered pf's
foreach g in xautoreduce(foreach f in p collect xreorder f)
  collect repartit g;


% ----------------------------------------------------------------------


symbolic procedure xautoreduce p;
% p is a list of reordered pf's, result is a list of reordered pf's
begin scalar xreducetime!*, % so reduction and auto-reduction times
                            % disjoint.
             q,s,x,t0;
 if !*xstats then xreducetime!* := t0 := time();
 while p do
   begin
%    s := xleast p; p := delete(s,p);
    s := car p; p := cdr p;     % seems to be faster
    x := xreduce(s,append(p,q));
    if x then q := x . q;
    if x neq s then
      <<p := append(p,q); q := nil>>;
   end;
 if !*xstats then xautoreducetime!* := xautoreducetime!* + time() - t0;
 return reverse q;
end;


% ----------------------------------------------------------------------


symbolic procedure xreduce0(s,p);
% s is pf, p is list of pf's, result is pf
% symbolic mode entry point to xreduce for unordered pf's

xreduce1(xreorder s,foreach f in p collect xreorder f);


% ----------------------------------------------------------------------


symbolic procedure xreduce(s,p);
% s is reordered pf, p is list of reordered pf's, result is reordered pf
begin scalar t0,xwedgestime!*; % so wedge and reduction times disjoint
 if !*xstats then t0 := time();
 s := xreorder xreduce1(s,p);
 if !*xstats then xreducetime!* := xreducetime!* + time() - t0;
 return s;
end;


% ----------------------------------------------------------------------


symbolic procedure xreduce1(s,p);
% s is reordered pf, p is list of reordered pf's, result is pf (NOT
% reordered)
if null s then nil else
begin scalar q,f,x,y,z;
 q := p;
 while s and q do
   begin
%    f := xleast q; q := delete(f,q);
    f := car q; q := cdr q;     % seems to be faster
    if x := xdivs(ldpffax s,ldpffax f) then
      <<y := partitsq!* simp list('wedge,x,prepsq !*pf2sq(lt f .+ nil));
        z := negsq quotsq(lc s,lc y);
        y := partitsq!* simp list('wedge,x,prepsq !*pf2sq f);
        % the repartit here needed because of ordering differences
        s := xreorder repartit addpf(s,multpfs(1 .* z .+ nil,y));
        q := p>>;
   end;
 return if s then addpf(lt s .+ nil,xreduce1(red s,p));
end;


% ----------------------------------------------------------------------


symbolic procedure xdivs(u,v);
% u,v are kernels, result is kernel
% seems to be faster than old one
begin scalar x;
 if setdiff(v,intersection(u,v)) then return nil;
 if (x := setdiff(u,v)) then return 'wedge . x
 else return 1;
end;


% ----------------------------------------------------------------------


symbolic procedure ldpffax u;
% returns list of factors in a wedge product (including single factors)

(if atom x or car x neq 'wedge then list x else cdr x) where x = ldpf u;


% ----------------------------------------------------------------------


symbolic procedure wedgepf(u,v);
% u and v are pf, result is pf
partitwedge list(mk!*sq !*pf2sq u,mk!*sq !*pf2sq v);


% ----------------------------------------------------------------------


symbolic procedure xreorder u;
% reorders pf u into (graded) lexicographical order. basic ordering
% of 1-forms determined by worderp.
sort(u,function xordpft);


% ----------------------------------------------------------------------


symbolic procedure xordpft(u,v);
% u and v are terms from a pf
apply(function xgradlexl,              % can put xlexl or xgradlexl here
      list(sort(ldpffax list u,function worderp),
           sort(ldpffax list v,function worderp)));


% ----------------------------------------------------------------------


symbolic procedure xlexl(u,v);
% u and v are sorted lists of factors from wedge products
if null v then t
else if null u then nil
else if car v = 1 then t
else if car u = 1 then nil
else if car u = car v then xlexl(cdr u,cdr v)
else worderp(car u,car v);


% ----------------------------------------------------------------------


symbolic procedure xgradlexl(u,v);
% u and v are sorted lists of factors from wedge products
if car v = 1 then t
else if car u = 1 then nil
else if length u = length v then xlexl(u,v)
else length u > length v;


% ----------------------------------------------------------------------


symbolic procedure xleast p;
% p is a list of reordered pf's, result is a reordered pf
begin scalar s,x,y;
 if p then s := car p;
 p := nil . p;
 while (p := cdr p) do
   if xordpft(lt s,lt car p) then s := car p;
 return s;
end;


% ----------------------------------------------------------------------


symbolic procedure xdegree qform;
% xdegree(qform: pform (prefix) )
% -> q: integer
%
% This procedure gives the degree of a homogeneous form.
% Can distinguish 0- and p-forms.
% Behaves erratically with inhomogeneous forms.

(if null x then 0 else x)
  where x = (deg!*form qf)
    where qf = if eqcar(qform,'!*sq) then prepsq cadr qform else qform;


% ----------------------------------------------------------------------


symbolic procedure xvarcheck(u,v);
% u is pf, v is list of kernels, result is list of kernels

if null u then v
else begin scalar vv;
 vv := setdiff(ldpffax u, 1 . v);
 v := append(v,vv);
 if fixp dimex!* and length v > dimex!* then
   rederr "more independent 1-forms than dimension of space in XIDEAL";
 foreach f in vv do if xdegree f > 1 then
   rederr "p-form with p > 1 in XIDEAL";
 return xvarcheck(red u,v);
end;


% ----------------------------------------------------------------------


symbolic procedure xhomogeneous u;
% u is pf, result is degree of u if homogeneous, otherwise nil

if null u then 0
else if length u = 1 then xdegree ldpf u
else (if d = xhomogeneous red u then d) where d = xdegree ldpf u;


% ----------------------------------------------------------------------


symbolic procedure xdegreecheck u;
% u is prefix
% result is non-nil if degree of u exceeds truncation in graded GB's

xtruncate!* and (xdegree u > xtruncate!*);


% ----------------------------------------------------------------------


xstatsinitialise(); % so that calls from other packages don't get
                    % non-numeric errors if xstats happens to be on.



% ----------------------------------------------------------------------


endmodule;


module xidealutils;

% Various list utilities

% ----------------------------------------------------------------------


symbolic procedure !*al2sl l;

% Converts algebraic mode lists to symbolic mode lists, with an error
% if conversion not possible

if atom l or car l neq 'list then typerr(l,'list) else
  for each j in cdr l collect if eqcar(j,'!*sq) then prepsq cadr j
                               else j;


% ----------------------------------------------------------------------


symbolic procedure deleteall(u,v);
% removes all top-level occurences of u from v
if null v then nil
else if car v = u then deleteall(u, cdr v)
     else car v . deleteall(u, cdr v);


% ----------------------------------------------------------------------


symbolic procedure purge l;
% removes extra copies of repeated elements from algebraic and symbolic
% lists
if null l then nil
else car l . purge deleteall(car l,cdr l);


% ----------------------------------------------------------------------


symbolic procedure strip0 l;
% removes all all NIL's and 0's from algebraic and symbolic lists
if null l then nil
else deleteall(0,deleteall(nil,l));


% ----------------------------------------------------------------------

endmodule;


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