reduce-historical/r38/packages/eds/systems.red

module systems;

% Operations on exterior differential systems

% Author: David Hartley

fluid '(kord!* xtruncate!* !*arbvars !*edssloppy cfrmcob!*);
global '(indxl!*);


symbolic procedure copyeds s;
   % s:eds -> copyeds:eds
   % Copy s to allow destructive operations using selectors
   append(s,{});

 

put('augment,'rtypefn,'getrtypecar);
put('augment,'edsfn,'augmenteds);

symbolic procedure augmenteds(s,u);
   % s:eds, u:prefix sys -> augmenteds:eds
   begin
   u := makelist getrlist u;
   u := !*a2sys u;
   s := augmentsys(s,u);
   foreach f in {'pfaffian,'closed,'quasilinear,'involutive} do
      rempropeds(s,f);
   return checkeds s; % removes all hidden flags, adds new rsx
   end;


symbolic procedure augmentsys(s,u);
   % s:eds, u:sys -> augmentsys:eds
   % Augment system by adding new forms, using ordering of s on input in
   % final sort. Doesn't change flags or check integrity.
   begin scalar c;
   s := copyeds s;
   eds_sys s := sortsys(union(u,eds_sys s),edscob s);
   return s;
   end;


put('quasilinear,'psopfn,'quasilineareval);

symbolic procedure quasilineareval s;
   % s:{eds} -> quasilineareval:0 or 1
   if edsp(s := reval car s) then
      if knowntrueeds(s,'quasilinear) or
      	 not knownfalseeds(s,'quasilinear) and
      	 edscall quasilinearp s then 1 else 0
   else typerr(s,'eds);


symbolic procedure quasilinearp s;
   % s:eds -> quasilinearp:bool
   % Test whether (closure of) system is quasilinear
   knowntrueeds(s,'quasilinear) or
      not knownfalseeds(s,'quasilinear) and
      if not normaledsp s then
      	 rerror(eds,000,{"System not in normal form in quasilinearp"})
      else if null scalarpart eds_sys s and
	      quasilinearsys(nonpfaffpart eds_sys closure s,prlkrns s)
       then << flagtrueeds(s,'quasilinear); t >>
      else
   	 << flagfalseeds(s,'quasilinear); nil >>;


symbolic procedure quasilinearsys(s,prl);
   % s:sys, prl:list of 1-form kernel -> quasilinearsys:bool
   % Systems with 0-forms are non-linear by definition here.
   null cadr lineargens(s,{},prl);
   

symbolic procedure lineargenerators s;
   % s:eds -> lineargenerators:eds
   % Makes linearly generated part of s explicitly linear.
   begin scalar p;
   p := pfaffpart eds_sys s;
   p := append(p,append(car q,cadr q)) 
      	   where q = lineargens(setdiff(eds_sys s,p),{},prlkrns s);
   if p = eds_sys s then return s;
   s := copyeds s;
   eds_sys s := p;
   return sorteds s;
   end;


symbolic procedure lineargens(s,g,prl);
   % s,g:sys, prl:list of 1-form kernel -> lineargens:{sys,sys}
   % g is a GB for a linear system, s is fully reduced wrt g. Returns as
   % linear a set of generators as possible.  For a linear system,
   % returns a linear set of generators.  Recursively checks if
   % non-linear part of s can be reduced mod linear part + g to give a
   % linear system.  Systems with 0-forms are non-linear by definition
   % here.
   begin scalar w,xtruncate!*; integer d;
   foreach f in s do
   << d := max(d,degreepf f);
      if degreepf f neq 0 and quasilinearpf(f,prl) then
 	 w := f . w >>;
   w := reversip w;
   s := setdiff(s,w);
   if null s then return {w,{}};
   if null w then return {{},s};
   xtruncate!* := d;
   g := xidealpf append(g,w);
   s := foreach f in s join
           if f := xreduce(f,g) then {f};
   return {append(w,car p),cadr p} where p = lineargens(s,g,prl);
   end;
   

symbolic procedure quasilinearpf(f,p);
   % f:pf, p:list of 1-form kernel -> quasilinearpf:bool
   % result is t if f is at most linear in p
   if null f then t
   else length intersection(wedgefax lpow f,p) <= 1
        and quasilinearpf(red f,p);


put('semilinear,'psopfn,'semilineareval);

symbolic procedure semilineareval s;
   % s:{eds} -> semilineareval:0 or 1
   if edsp(s := reval car s) then
      if edscall semilinearp s then 1 else 0
   else typerr(s,'eds);


symbolic procedure semilinearp s;
   % s:eds -> semilinearp:bool
   % Test whether (closure of) system is semilinear
   if not normaledsp s then nil
   else if !*edssloppy then edscall quasilinearp s
   else semilinearsys(nonpfaffpart eds_sys edscall closure s,prlkrns s);


symbolic procedure semilinearsys(s,prl);
   % s:sys, prl:list of 1-form kernel -> semilinearsys:bool
   % 0-forms are non-linear by definition here.
   null s or
      degreepf car s neq 0 and
      semilinearpf(car s,prl) and
      semilinearsys(cdr s,prl);

   
symbolic procedure semilinearpf(f,p);
   % f:pf, p:list of 1-form kernel -> semilinearpf:bool
   % Works when xvars!* allows 0-forms as well - used in solvegraded.
   % result is t if f is at most linear in p with constant coefficient
   null f or
      (l = 0 or
       	 l = 1 and
 	 cfrmconstant numr lc f and
 	 cfrmconstant denr lc f
	    where l = length foreach k in wedgefax lpow f join
 	       	     	      	if k memq p then {k}) and
      semilinearpf(red f,p);


put('pfaffian,'psopfn,'pfaffianeval);

symbolic procedure pfaffianeval s;
   % s:{eds} -> pfaffianeval:0 or 1
   if edsp(s := reval car s) then
      if knowntrueeds(s,'pfaffian) or
      	 not knownfalseeds(s,'pfaffian) and
      	 edscall pfaffian s then 1 else 0
   else typerr(s,'eds);


symbolic procedure pfaffian s;
   % s:eds -> pfaffian:bool
   knowntrueeds(s,'pfaffian) or
      not knownfalseeds(s,'pfaffian) and
      if not normaledsp s then
      	 rerror(eds,000,{"System not in normal form in pfaffian"})
      else if pfaffsys eds_sys s then
      << flagtrueeds(s,'pfaffian); t>>
      else
      << flagfalseeds(s,'pfaffian); nil>>;


symbolic procedure pfaffsys s;
   % s:sys -> pfaffsys:bool
   % Systems with 0-forms are non-Pfaffian by definition here.
   begin scalar p,xtruncate!*; integer d;
   if scalarpart s then return nil;
   foreach f in s do
   << d := max(d,degreepf f);
      if degreepf f = 1 then p := f . p >>;
   s := setdiff(s,p);
   if null s then return t;
   if null p then return nil;
   xtruncate!* := d;
   p := xidealpf foreach f in p collect xreduce(exdfpf f,p);
   while s and null xreduce(car s,p) do s := cdr s;
   return null s;
   end;


put('closure,'edsfn,'closure);
put('closure,'rtypefn,'getrtypecar);

symbolic procedure closure s;
   % s:eds -> closure:eds
   begin scalar p,sys,s0; integer d;
   if knowntrueeds(s,'closed) then return s;
   if s0 := geteds(s,'closure) then return s0;
%%%   if not normaledsp s then
%%%      rerror(eds,000,{"System not in normal form in closure"});
%%%   if scalarpart eds_sys s then
%%%      rerror(eds,000,{"Closure with 0-forms not yet implemented"});
   if scalarpart eds_sys s then
      lprim {"0-forms in closure: result may not be closed"};
   d := length eds_ind s;
   p := solvedpart eds_sys s;
   sys := foreach f in eds_sys s join
             if degreepf f < d and 
		(f := xreduce(xreorder exdfpf f,p)) then {f};
   if null sys then return <<flagtrueeds(s,'closed); s>>;
   s0 := augmentsys(s,sys);
   if pfaffpart sys then rempropeds(s0,'solved);
   flagtrueeds(s0,'closed);
   s0 := normaleds s0;        % might add 0-forms or become inconsistent
   return if emptyedsp s0 then s0
      else if scalarpart eds_sys s0 then s0
      else <<puteds(s,'closure,s0); s0>>;
   end;

flag('(closure),'hidden);


% symbolic operator closed;
% symbolic procedure closed u;
%    % u:eds|rlist of prefix|prefix -> closed:bool
%    % True if u is a closed eds, a closed system of forms or a closed
%    % form
%    if edsp u then
%       knowntrueeds(u,'closed) or edscall closededs u
%    else if rlistp u then
%       closedsys foreach f in getrlist u collect xpartitop f
%    else null exdfpf xpartitop u;

% flag('(closed),'boolean);

% symbolic procedure closededs s;
%    % s:eds -> closededs:bool
%    knowntrueeds(s,'closed) or
%    if closedsys eds_sys s then
%    << flagtrueeds(s,'closed); t>>;



put('closed,'psopfn,'closedeval);

symbolic procedure closedeval s;
   % s:{eds} -> closedeval:0 or 1
   if edsp(s := reval car s) then
      if knowntrueeds(s,'closed) or
      	 not knownfalseeds(s,'closed) and
         edscall closed s then 1 else 0
   else if rlistp s then
      if closedsys foreach f in getrlist s collect xpartitop f then 1
      else 0
   else if null exdfpf xpartitop s then 1 else 0;


symbolic procedure closed s;
   % s:eds -> closed:bool
   knowntrueeds(s,'closed) or
      not knownfalseeds(s,'closed) and
      if closedsys eds_sys s then
      << flagtrueeds(s,'closed); t>>
      else
      << flagfalseeds(s,'closed); nil>>;


symbolic procedure closedsys s;
   % s:sys -> closedsys:bool
   begin scalar p,xtruncate!*; integer d;
   foreach f in s do
   << d := max(d,1 + degreepf f);
      f := xreduce(exdfpf f,s);
      if f then p := f . p >>;
   if null p then return t;
   xtruncate!* := d;
   s := xidealpf s;
   while p and null xreduce(car p,s) do p := cdr p;
   return null p;
   end;


symbolic operator frobenius;
symbolic procedure frobenius u;
   % u:eds|rlist of prefix|prefix -> frobenius:bool
   % True if u is an eds or list of forms generated by 1-forms
   % satisfying the Frobenius test
   if edsp u then
      null nonpfaffpart eds_sys u and
      null nonpfaffpart eds_sys edscall closure u
   else if rlistp u then
      frobeniussys foreach f in getrlist u collect xpartitop f
   else rerror(eds,000,"Invalid argument to frobenius");

flag('(frobenius),'boolean);


symbolic procedure frobeniussys s;
   % s:sys -> frobeniussys:bool
   begin scalar p;
   p := pfaffpart s;
   s := union(foreach f in p collect exdfpf f,setdiff(s,p));
   p := xautoreduce p;
   while s and null xreduce(car s,p) do s := cdr s;
   return null s;
   end;


put('cauchy_system,'rtypefn,'quotelist);
put('cauchy_system,'listfn,'evalcauchysys);

symbolic procedure evalcauchysys(u,v);
   % u:{prefix}, v:bool -> evalcauchysys:rlist
   if xedsp(u := reval car u) then
      evalcartansys({edscall closure u},v)
   else if rlistp u then
      evalcartansys({append(u,foreach f in cdr u collect aeval{'d,f})},v)
   else
      evalcartansys({makelist {u,aeval{'d,u}}},v);


put('cartan_system,'rtypefn,'quotelist);
put('cartan_system,'listfn,'evalcartansys);

symbolic procedure evalcartansys(u,v);
   % u:{prefix}, v:bool -> evalcartansys:rlist
   if xedsp(u := reval car u) then
      if edsp u then !*sys2a1(edscall cartansyseds u,v)
      else makelist for each s in cdr u
		       collect !*sys2a1(edscall cartansyseds u,v)
   else if rlistp u then
      !*sys2a1(cartansys !*a2sys u,v)
   else
      !*sys2a1(cartansyspf xpartitop u,v);


symbolic procedure cartansys u;
   % u:list of pf -> cartansys:list of pf
   begin scalar xtruncate!*;
   xtruncate!* := eval('max.foreach f in u collect degreepf f);
   xtruncate!* := xtruncate!* - 1;
   u := xidealpf u;
   return reversip xautoreduce purge foreach f in u join cartansyspf f;
   end;


symbolic procedure cartansyspf f;
   % f:pf -> cartansyspf:list of pf
   begin scalar x,p,q,z;
   if degreepf f = 1 then return {f};
   while f do
      begin
      p := wedgefax lpow f;
      foreach k in p do
         if not((q := delete(k,p)) member z) then
	 << z := q . z;
	    x := xcoeff(f,q) . x >>;
      f := red f;
      end;
   return reverse xautoreduce purge x;
   end;


symbolic procedure cartansyseds s;
   % s:eds -> cartansyseds:sys
   cartansys eds_sys s;


put('linearise,'edsfn,'lineariseeds);
put('linearise,'rtypefn,'quoteeds);
put('linearize,'edsfn,'lineariseeds);
put('linearize,'rtypefn,'quoteeds);
flag('(linearise linearize),'nospread);

symbolic procedure lineariseeds u;
   % u:{eds[,rlist]} -> lineariseeds:eds
   begin scalar x;
   if null u or length u > 2 then
      rerror(eds,000,{"Wrong number of arguments to linearise"});
   if cdr u then x := !*a2sys cadr u;
   if nonpfaffpart x then typerr(cadr u,"integral element");
   return edscall linearise(car u,x);
   end;


symbolic procedure linearise(s,x);
   % s:eds, x:sys -> linearise:eds
   % x is an integral element of s, result is linearisation of s at x
   % in original cobasis.
   if quasilinearp s then lineargenerators s else
   begin scalar xx,q,prl;
   s := copyeds closure s;
   x := xreordersys x;
   q := nonpfaffpart eds_sys s;
   prl := prlkrns s;
   % pick out those products which occur
   xx := purge foreach f in q join
      	    foreach k in xpows f join
	       nonlinfax intersection(wedgefax k,prl);
   % form the relevant poducts from x
   x := pair(lpows x,x);
   xx := foreach pr in xx collect 
      	    wedgepf(cdr atsoc(car pr,x),cdr atsoc(cadr pr,x));
   % reduce the system mod x^x
   eds_sys s := append(setdiff(eds_sys s,q),
      	        foreach f in q join if f := xreduce(f,xx) then {f});
   flagtrueeds(s,'quasilinear);
   return s;
   end;


symbolic procedure nonlinfax l;
   % l:list of kernel -> nonlinfax:list of list of 2 kernel
   % Collect elements of l pairwise, discarding any odd element.
   if length l > 1 then {car l,cadr l} . nonlinfax cddr l;

%% symbolic procedure linearise(s,x);
%%    % s:eds, x:sys -> linearise:eds
%%    % x is an integral element of s, result is linearisation of s at x
%%    % in original cobasis.
%%    % NB Changes background coframing.
%%    if quasilinearp s then lineargenerators s else
%%    begin scalar s1;
%%    s1 := linearise0(s,x);
%%    x := cadr s1;
%%    s1 := car s1;
%%    return xformeds0(s1,x,setdiff(edscob s,edscob s1));
%%    end;
%% 
%% 
%% symbolic procedure linearise0(s,x);
%%    % s:eds, x:sys -> linearise0:{eds,xform}
%%    % x is an integral element of s, result is linearisation of s at x
%%    % in a cobasis based on x, together with transform required for
%%    % original cobasis. The structure equations are NOT updated.
%%    % NB Changes background coframing.
%%    begin scalar c,y,prl;
%%    c := foreach f in x collect mkform!*(intern gensym(),1);
%%    x := pair(c,x);
%%    y := invxform x;
%%    s := copyeds closure s;
%%    s := tmpind xformeds0(s,y,c);
%%    x := append(x,cadr s);
%%    s := car s;
%%    prl := prlkrns s;
%%    eds_sys s := foreach f in eds_sys s join
%%       	       	   if degreepf f < 2 then {f}
%% 		   else if inhomogeneouspart(f,prl) then
%% 		      typerr(!*sys2a foreach p in x collect cdr p,
%% 			     "integral element")
%% 		   else if f := linearpart(f,prl) then {f};
%%    flagtrueeds(s,'quasilinear);
%%    return {s,x};
%%    end;


put('one_forms,'rtypefn,'quotelist);
put('one_forms,'listfn,'oneformseval);

symbolic procedure oneformseval(u,v);
   % u:{xeds|rlist}, v:bool -> oneformseds:rlist
   if edsp(u := reval car u) then
      !*sys2a1(pfaffpart eds_sys u,v)
   else if xedsp u then
      makelist foreach s in getrlist u collect
	 !*sys2a1(pfaffpart eds_sys s,v)
   else
      makelist foreach f in getrlist u join
	 if reval{'exdegree,f}=1 then {f};


put('zero_forms,'rtypefn,'quotelist);
put('zero_forms,'listfn,'zeroformseval);
put('nought_forms,'rtypefn,'quotelist);
put('nought_forms,'listfn,'zeroformseval);

symbolic procedure zeroformseval(u,v);
   % u:{xeds|rlist}, v:bool -> zeroformseval:rlist
   if edsp(u := reval car u) then
      !*sys2a1(scalarpart eds_sys u,v)
   else if xedsp u then
      makelist foreach s in getrlist u collect
	 !*sys2a1(scalarpart eds_sys s,v)
   else
      makelist foreach f in getrlist u join
	 if reval{'exdegree,f}=0 then {f};

endmodule;

end;