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

module edsuser;

% Miscellaneous user functions

% Author: David Hartley


fluid '(alglist!* subfg!* !*arbvars xvars!*);

global '(indxl!*);


put('index_expand,'rtypefn,'quotelist);
put('index_expand,'listfn,'indexexpandeval0);

symbolic procedure indexexpandeval0(u,v);
   % u:list of prefix, v:bool -> indexexpandeval0:prefix list
   % kludge to add v argument to index_expand listfn
   makelist foreach p in getrlist indexexpandeval u collect reval1(p,v);


symbolic procedure indexexpandeval u;
   % u:list of prefix -> indexexpandeval:prefix list
   if length u neq 1 then 
      rerror(eds,000,"Wrong number of arguments to index_expand")
   else if rlistp(u := reval car u) then
      makelist purge foreach x in cdr u join cdr indexexpandeval {x}
   else makelist indexexpand u;


symbolic procedure indexexpand u;
   % u:prefix -> indexexpand:list of prefix
   % u has always been reval'd, so there is no need to expand
   % summations.
   if eqexpr u then indexexpandeqn u else
   if boolexpr u then indexexpandbool u else
   begin scalar i,v,alglist!*;
   u := simp!* u;
   % Expand free indices (put them into some order for safety)
   i := idsort purge if numr u and not domainp numr u then
      	       	 flatindxl allind !*t2f lt numr u;
   v := foreach j in mkaindxc(i,nil) join
      	   if numr(j := subfreeindices(numr u,pair(i,j))) then
	      {absf numr j ./ denr j};  % nprimitive too?
   return for each q in purge v collect mk!*sq multsq(q,1 ./ denr u);
   end;


symbolic procedure indexexpandeqn u;
   % u:rule|equation -> indexexpandeqn:list of rule|equation
   begin scalar i,v,lh,rh;
      	 scalar alglist!*;
   % Expand free indices on lhs (put them into some order for safety)
   lh := reval cadr u where subfg!* = nil; % avoid let rules
   i := idsort purge flatindxl allindk lh;
   rh := aeval caddr u;
   v := foreach j in mkaindxc(i,nil) join
      if j := subfreeindeqn({car u,lh,rh},pair(i,j)) then {j};
   % Remove duplicates
   i := {};
   v := foreach r in v join
      if not(cadr r member i) then
      	 << i := cadr r . i; {r} >>;
   return v;
   end;


symbolic procedure subfreeindeqn(u,l);
   % u:rule|equation, l:alist -> subfreeindeqn:rule|equation|nil
   % Make index substitution l in u. Only index symmetry simplifications
   % are allowed, so the lhs can either vanish (nil returned) or acquire
   % an overall sign (sign transferred to rhs).
   begin scalar lh,rh;
   lh := subfreeindk(cadr u,l);
   if null atomf lh then lh := revop1 lh; % gets done in rule!-list
					  % anyway
   lh := reval lh where subfg!* = nil; % avoid let rules;
   if lh = 0 then return nil;
   rh := simp!* caddr u;
   rh := quotsq(subfreeindices(numr rh,l),subfreeindices(denr rh,l));
   if eqcar(lh,'minus) then
   << lh := cadr lh;
      rh := negsq rh >>;
   return {car u,lh,mk!*sq rh};
   end;


symbolic procedure boolexpr u;
   % u:any -> boolexpr:bool
   not atom u and flagp(car u,'boolean);


symbolic procedure indexexpandbool u;
   % u:prefix -> indexexpandbool:list of prefix
   begin scalar i,v,alglist!*;
   % Expand free indices on lhs (put them into some order for safety)
   i := idsort purge flatindxl allindk u;
   v := foreach j in mkaindxc(i,nil) collect
	   car u .
	      foreach a in cdr u collect reval subfreeindk(a,pair(i,j));
   return purge v;
   end;


symbolic procedure subfreeindices(u,l);
   % u:sf, l:a-list -> subfreeindices:sq
   % Discriminates indices from variables, modified from EXCALC's
   % subfindices to go inside operators other than EXCALC's.
   begin scalar alglist!*;
   return 
      if domainp u then u ./ 1
      else 
	 addsq(
	    multsq(if atom mvar u then
	       	      !*p2q lpow u
	       	   else if sfp mvar u then
		      exptsq(subfreeindices(mvar u,l),ldeg u)
		   else 
		      exptsq(simp subfreeindk(mvar u,l),ldeg u),
		   subfreeindices(lc u,l)),
	    subfreeindices(red u,l))
   end;


symbolic procedure subfreeindk(u,l);
   % u:kernel, l:a-list -> subfreeindk:kernel
   % Extends subindk to indexed variables
   if atom u then u
   else if flagp(car u,'indexvar) then
      car u . subla(l,cdr u)
   else subindk(l,u);


put('linear_divisors,'rtypefn,'quotelist);
put('linear_divisors,'listfn,'lineardivisors);

symbolic procedure lineardivisors(u,v);
   % u:{prefix}, v:bool -> lineardivisors:prefix list
   makelist foreach f in lineardivisorspf xpartitop reval car u collect
      !*pf2a1(f,v);


symbolic procedure lineardivisorspf f;
   % f:pf -> lineardivisorspf:list of pf
   begin scalar x,g,v;
   foreach p in xpows f do x := union(wedgefax p,x);
   foreach k in x do
   << v := intern gensym() . v;
      g := addpf(k .* !*k2q car v .+ nil,g)>>;
   x := edssolve(xcoeffs wedgepf(g,f),v);
   if length x neq 1 or caar x neq t then
      errdhh "Bad solve result in lineardivisorspf";
   x := cadar x;
   v := updkordl v;
   g := numr subf(numr !*pf2sq g,x);
   x := {};
   while g do
   << x := xpartitsq(lc g ./ 1) . x; 
      g := red g >>;
   setkorder v;
   return reverse xautoreduce x;
   end;
   

symbolic procedure xdecomposepf f;
   % f:pf -> xdecomposepf:list of pf
   begin scalar x;
   x := lineardivisorspf f;
   if length x = degreepf f then return reverse x;
   end;


put('exfactors,'rtypefn,'quotelist);
put('exfactors,'listfn,'exfactors);

symbolic procedure exfactors(u,v);
   % u:{prefix}, v:bool -> exfactors:prefix list
   makelist foreach f in xfactorspf xpartitop reval car u collect
      !*pf2a1(f,v);


symbolic procedure xfactorspf f;
   % f:pf -> xfactorspf:list of pf
   begin scalar x;
   x := lineardivisorspf f;
   f := xreduce(f,foreach g in x collect 
      	       	     addpf(1 .* (-1 ./ 1) .+ nil,g));
   return
      if f = !*k2pf 1 then reverse x
      else f . reverse x;
   end;


symbolic operator exact;
symbolic procedure exact u;
   % u:prefix -> exact:bool
   % True if u is an exact pform kernel
   eqcar(u,'d);

flag('(exact),'boolean);


put('derived_system,'rtypefn,'getrtypecar);
put('derived_system,'edsfn,'deriveeds);
put('derived_system,'listfn,'derivelist);

symbolic procedure derivelist(u,v);
   % u:{xeds|rlist}, v:bool -> derivelist:rlist
   !*sys2a1(derive !*a2sys reval car u,v) where xvars!* = nil;


symbolic procedure deriveeds s;
   % s:eds -> deriveeds:eds
   begin
   s := copyeds s;
   if pfaffian s then
      eds_sys s := derive pfaffpart eds_sys s
   else rerror(eds,000,"non-Pfaffian system in derived_system");
   return s;
   end;


symbolic procedure derive s;
   % s:sys -> derive:sys
   begin scalar c,f;
   if null s then s;
   s := xautoreduce s;
   c := for each f in s collect 
      	 if degreepf f = 1 then intern gensym()
	 else rerror(eds,000,"non-Pfaffian system in derived_system");
   f := zippf(s,foreach k in c collect !*k2q k);
   s := edssolve(xcoeffs xreduce(exdfpf f,s),c);
   if length s neq 1 or null caar s then
      errdhh{"Bad solve result in derive:",s};
   s := cadr car s;
   f := pullbackpf(f,s);
   c := setdiff(c,foreach m in s collect car m);
   f := xrepartit f where xvars!* = c;
   return for each x in reverse c collect xcoeff(f,{x});
   end;


symbolic procedure allcoords f;
   % f:prefix -> allcoords:list of kernel
   % Pick up 0-form kernels in f
   makelist purge 
      foreach k in (xvarspf xpartitop f where xvars!* = t) join
      	 if xdegree k = 0 and 
	    not assoc(k,depl!*) and
 	    not eqcar(k,'partdf) then {k}
	 else if (xdegree k = 1) and exact k then {cadr k};

endmodule;

end;