Check-ins Using Download Hex

Artifact 6a14be384e50bad0222d28a6f0af932c81e7a806910bb5cd35737d03b6b380c0:

module edssolve;

% Specialised solvers for EDS

% Author: David Hartley

Comment. The EDS solve routines are interfaces to the general REDUCE
solver, presenting the results in a form more useful in EDS. The most
primitive is edssolve, which basically just turns off arbvars and
reorganises the answer, and the most sophisticated is edssolvegraded,
which tries to solve the principal part first if the equations are
semilinear. This is useful since the principal part is often less
complex than the rest. Solutions are returned as rmaps, giving
information about both the solution and its applicability.

endcomment;


fluid '(!*edsverbose !*edsdebug !*arbvars !*varopt !*groebopt
        !*solveinconsistent depl!* cfrmcrd!* cfrmrsx!* xvars!*
	!*edssloppy !*edsdisjoint);


symbolic procedure edssolvegraded(xl,vl,rsx);
   % xl:list of sq, vl:list of list of kernel, rsx:rsx
   % -> edssolvegraded:list of (t.rmap)|(nil.list of sq)
   (begin scalar fl,z;
   z := foreach l in vl join append(l,{});
   if vl and semilinearsql(xl,car vl,if !*edssloppy then car vl else z)
     then return edssolvesemi(xl,vl,rsx);
   foreach l on vl do
      foreach y in car l do foreach ll in cdr l do
 	 foreach x in ll do depend1(y,x,t);
   xl := edssolve(xl,z);
   fl := {};
   xl := foreach s in xl join
      if null car s  then <<fl := s . fl; nil>>
      else if not member(0,z := pullbackrsx(rsx,cadr s)) then
 	 {{cadr s,purgersx append(z,caddr s)}};
   return append(reverse fl,foreach s in edsdisjoin xl collect t . s);
   end) where depl!* = depl!* ;


symbolic procedure semilinearsql(xl,vl,kl);
   % xl:list of sq, vl,kl:list of kernel -> semilinearsql:bool
   null xl
      or semilinearsq(car xl,vl,kl) and semilinearsql(cdr xl,vl,kl);

   
symbolic procedure semilinearsq(x,vl,kl);
   % x:sq, vl,kl:list of kernel -> semilinearsq:bool
   % True if x is linear in vl with coefficients independent of kl.
   % Assumes that kl includes vl.
   semilinearpf0(xpartitsq x,vl,kl) where xvars!* = vl;


symbolic procedure semilinearpf0(f,vl,kl);
   % f:pf, vl,kl:list of kernel  -> semilinearpf0:bool
   % Works when xvars!* allows 0-forms as well - used in solvegraded.
   % True if x is linear in vl with coefficients independent of kl.
   % Assumes that kl includes vl.
   null f or
      (lpow f = 1 and 
       freeoffl(denr lc f,vl) and   % because vl might occur inside
				    % operators
       freeoffl(numr lc f,vl)
       	  or
       length wedgefax lpow f = 1 and
       freeoffl(denr lc f,kl) and   % because vl might occur inside
				    % operators
       freeoffl(numr lc f,kl)) and
      semilinearpf0(red f,vl,kl);


symbolic procedure edssolvesemi(xl,vl,rsx);
   % xl:list of sq, vl:list of list of kernel, rsx:rsx
   % -> edssolvesemi:list of (t.rmap)|(nil.list of sq)
   % xl is known to be semilinear
   begin scalar sl,nl;
   xl := edspartsolve(xl,car vl);
   foreach x in car xl do
      if not(car x memq cadr xl) then
	 sl := (car x . mk!*sq subs2 subsq(simp!* cdr x,caddr xl)) . sl
      else if numr
      << x := addsq(negsq !*k2q car x,simp!* cdr x);
	 x := subs2 subsq(x,caddr xl) >>
      then nl := x . nl;
   sl := !*map2rmap reversip sl;
   if 0 member (rsx := pullbackrsx(rsx,car sl)) then return nil;
   rsx := purgersx append(rsx,cadr sl); sl := car sl;
   if null nl then return {t . {sl,rsx}};
   xl := edssolvegraded(nl,cdr vl,rsx);
   return foreach s in xl collect
      if null car s then
	 nil . append(cdr s,foreach x in sl collect 
  	       	     	       addsq(negsq !*k2q car x,simp!* cdr x))
      else
	 t . {append(pullbackmap(sl,cadr s),cadr s),caddr s};
   end;


symbolic procedure edssolve(xl,vl);
   % xl:list of sq, vl:list of kernel -> edssolve: list of tag.value
   % where tag.value is one of
   % 	t.rmap		an explicit solution
   % 	nil.list of sq	a (partial) system which couldn't be solved
   % This is an interface to the REDUCE solve routines, for use by other
   % eds routines. It switches off arbvars and returns the result in a
   % form more easily used in eds. If the system is inconsistent, an
   % empty list {} is returned.
   begin scalar kl,sl,msg;
         scalar !*arbvars;  % stop arbcomplex's being generated
   msg := {"Solving",length xl,"equations in",length vl,"variables"};
   foreach q in xl do kl := union(kernels numr q,kl);
   kl := expanddependence kl;
   vl := intersection(vl,kl);	% must preserve order of vl
   msg := append(msg,{{length vl,"present"}}); % {a,b} prints as (a b)
   edsdebug(msg,nil,nil);
   sl := edssolvesys(foreach q in xl collect numr q,vl);
   if eqcar(sl,'inconsistent) then
      sl := {}
   else if eqcar(sl,'failed) or eqcar(sl,nil) then
      sl := {nil . xl}
   else if eqcar(sl,t) then
      sl := foreach s in cdr sl join
      	 if not explicitsolutionp s then
            {nil . foreach pr in pair(cadr s,car s) collect
               addsq(negsq !*k2q car pr,cdr pr)}
      	 else % need to reorder return value from edssolvesys!
	    {t . !*map2rmap pair(cadr s,for each q in car s
					   collect mk!*sq reordsq q)}
   else errdhh{"Unexpected result from solvesys:",sl};
   return sl;
   end;


symbolic procedure expanddependence kl;
   % kl: list of kernel -> expanddependence:list of kernel
   % expand kl recursively to include all kernels explicitly appearing
   % as operator arguments. Should we check implicit dependence also?
   if null kl then nil
   else if atomf car kl then
      car kl . expanddependence cdr kl
   else % must be operator, so add all arguments to list
      car kl . expanddependence append(cdar kl,cdr kl);


symbolic procedure edssolvesys(xl,vl);
   % xl:list of sf, vl:list of kernel -> edssolvesys: list of tag.value
   % where tag.value is given in solve.red. This just calls solvesys.
   solvesys(xl,vl);


symbolic procedure explicitsolutionp s;
   % s:solve solution -> explicitsolutionp:bool
   not smember('root_of,s) and not smember('one_of,s);


symbolic procedure edspartsolve(xl,vl);
   % x:list of sq,
   % vl:list of kernel -> edspartsolve:{map,list of kernel,map}
   % Solves the equations xl for the variables vl by splitting into
   % homogeneous and inhomogeneous parts.  The results are only
   % guaranteed for linear equations whose coefficient are independent
   % of the inhomogeneous parts.  The solution is returned in terms of
   % some temporary variables representing the inhomogeneous parts.
   % The temporary variables are given in the original variables by the
   % third return value.
   (begin scalar al,il,l;
%	 scalar depl!*;	    % preserve dependencies
   xl := splitlinearequations(xl,vl);
   xl := foreach p in xl collect
      	    if null numr cdr p then
	       car p
	    else if (l := mk!*sq cdr p member il) then
	       addsq(car p,!*k2q nth(al,1 + length il - length l)) 
	    else
            << al := intern gensym() . al;
	       il := mk!*sq cdr p . il;
	       addsq(car p,!*k2q car al) >>;
   al := reversip al; il := reversip il;
   foreach y in vl do
      foreach z in al do depend1(y,z,t);
   xl := edssolve(xl,append(vl,al));
   if length xl neq 1 or null car(xl := car xl) then
      errdhh{"Bad solution in edspartsolve",xl};
   return {cadr xl,al,pair(al,il)};
   end) where depl!* = depl!*;	    % preserve dependencies


symbolic procedure splitlinearequations(v,c);
   % v:list of sq, c:list of kernel
   % -> splitlinearequations:{list of sq.sq}
   % Splits rational expressions in v into homogeneous and
   % inhomogeneous parts wrt variables in c.
   begin scalar ok,g,h;
         scalar !*exp;
   !*exp := t;
   ok := setkorder c;
   v := foreach q in v collect
        << g := h := nil;
      	   while not domainp f and mvar f memq c do
      	      << g := lt f . g; f := red f >>;
	   foreach u in g do h := u .+ h;
      	   cancel(h ./ d) . cancel(f ./ d) >>
	 where f = reorder numr q,
	       d = reorder denr q;
   setkorder ok;
   return foreach p in v collect reordsq car p . reordsq cdr p;
   end;


symbolic procedure edsgradecoords(crd,jet0);
   % crd,jet0:list of kernel -> edsgradecoords:list of list of kernel
   % grade coordinates according to jet order, highest jet first
   if null jet0 then {crd}
   else
      begin scalar u;
      foreach c in crd do
	 begin scalar j0,c0;
	 j0 := jet0;
	 while j0 and null(c0 := splitoffindices(car j0,c)) do
	    j0 := cdr j0;
	 if j0 := assoc(length c0,u) then nconc(j0,{c})
	 else u := (length c0 . {c}) . u;
	 end;
      u := sort(u,function(lambda x,y; car x > car y));
      return foreach v in u collect cdr v;
      end;



Comment. The routine solvepfsys tries to bring a system into solved form
in the current environment specified by cfrmcrd!* and cfrmrsx!*.

endcomment;

symbolic procedure solvepfsys sys;
   % sys:sys -> solvepfsys:{sys,sys}
   % Bring sys into solved form as far as possible.
   solvepfsys1(sys,{});

symbolic procedure solvepfsys1(sys,vars);
   % sys:sys, vars:list of lpow pf -> solvepfsys1:{sys,sys}
   % Solve sys for vars. If vars is {} then solve for anything. Kernel
   % ordering changed so that solved set comes ahead of rest. Ordering
   % within solved set and others is unchanged. The solved part is
   % returned sorted in decreasing order of lpow.  NBB. Kernel ordering
   % changed!!
   begin scalar ok,sl,nl;
   % save old kernel ordering
   ok := updkordl nil;
   nl := foreach f in sys collect subs2pf!* f;
   % First try for constant coefficients
   if nl then
      begin
%      nl := solvepfsyseasy(nl,vars,'domainp);
      nl := solvepfsyseasy(nl,vars,'cfrmconstant);
      if car nl then
 	 edsdebug("Solved with constant coefficients:",car nl,'sys);
      if cadr nl then
 	 edsdebug("Unsolved with constant coefficients:",cadr nl,'sys);
      sl := append(sl,car nl);
      nl := cadr nl;
      end;
   % If that's not enough, try for nowhere-zero coefficients
   if nl then
      begin
      nl := solvepfsyseasy(nl,vars,'cfrmnowherezero);
      if car nl then
 	 edsdebug("Solved with nonzero coefficients:",car nl,'sys);
      if cadr nl then
 	 edsdebug("Unsolved with nonzero coefficients:",cadr nl,'sys);
      sl := append(sl,car nl);
      nl := cadr nl;
      end;
   % If that's not enough, try Cramer's rule.
   if nl then
      begin
      nl := solvepfsyshard(nl,vars);
      if car nl then
 	 edsdebug("Solved with Cramer's rule:",car nl,'sys);
      if cadr nl then
 	 edsdebug("Unsolved with Cramer's rule:",cadr nl,'sys);
      sl := append(sl,car nl);
      nl := cadr nl;
      end;
   % Fix up kernel ordering and back-substitute solved forms
   setkorder ok;
   updkordl lpows sl;
   sl := xautoreduce1 xreordersys sl;
   nl := xreordersys nl;
   % Block structure may have messed up order of solved kernels
   sl := sortsys(sl,ok);
   updkordl lpows sl;
   return {sl,nl};
   end;


symbolic procedure solvepfsyseasy(sys,vars,test);
   % sys:sys, vars:list of lpow pf,
   % test:function -> solvepfsyseasy:{sys,sys}
   % Recursively bring sys into weakly solved form as far as possible
   % just by changing kernel ordering and normalising.  The solved part
   % is in normalised upper triangular form, and korder is completely
   % messed up - results must be reordered under lpows solved part.
   begin scalar sl,nl,kl;
   kl := purge foreach f in sys join xsolveables(f,test);
   kl := sort(if vars then intersection(kl,vars) else kl,'termordp);
   if null kl then return {{},sys};
   updkordl kl;
   foreach f in sys do
      if (f := xreduce(xreorder f,sl)) and apply1(test,numr lc f) then
         sl := xnormalise f . sl
      else if f then
	 nl := f . nl;
   sl := reversip sl; % sl in upper triangular form
   nl := reversip foreach f in nl join 
      if f := subs2pf!* xreduce(f,sl) then {f};
   if null sl or null nl then return {sl,nl};
   nl := solvepfsyseasy(nl,vars,test);
   return {append(sl,car nl),cadr nl};
   end;


symbolic procedure xsolveables(f,test);
   % f:pf, test:function -> xsolveables:list of lpow pf
   % All powers in f whose coefficient numerators satisfy test
   if null f then nil
   else if apply1(test,numr lc f) then lpow f . xsolveables(red f,test)
   else xsolveables(red f,test);


symbolic procedure solvepfsyshard(sys,vars);
   % sys:sys, vars:list of lpow pf -> solvepfsyshard:{sys,sys}
   % Bring sys into weakly solved form by Cramer's rule.  The solved
   % part is in normalised upper triangular form, and korder is
   % completely messed up - results must be reordered under lpows
   % solved part.  Allowing !*edssloppy to work here for the first time
   % minimises the number of restrictions added.
   begin scalar sl,nl,kl,w;
   if null sys then return {{},{}};
   if vars then updkordl vars;
   w := !*k2pf 1;
   foreach f in sys do
      if f := subs2pf wedgepf(if vars then xreorder f else f,w) then 
	 if null vars or subsetp(wedgefax lpow f,vars) then
	    w := f
	 else
	 << if degreepf f neq 1 then
	       f := xcoeff(f,intersection(wedgefax lpow f,vars));
	    nl := multpfsq(f,invsq lc w) . nl >>; % exact divisions
   kl := xsolveables(w,'cfrmnowherezero);
   if null kl and !*edssloppy then kl := xpows w;
   if vars then while kl and not subsetp(wedgefax car kl,vars) do
	 kl := cdr kl;
   if null kl or (car kl = 1) then return {{},sys};
   kl := wedgefax car kl;
   if !*edssloppy then
      cfrmrsx!* := xpartitsq(numr lc xcoeff(w,kl) ./ 1) . cfrmrsx!*
	       	     	 where xvars!* = cfrmcrd!*;
   sl := foreach k in kl collect xcoeff(w,delete(k,kl));
   updkordl kl;
   return {foreach f in sl collect xnormalise xreorder f,
      	   foreach f in nl collect xrepartit!* f};
   end;

endmodule;

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