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

module edspatch;

% Various patches for other parts of Reduce.

% Author: David Hartley

fluid '(!*edsverbose !*edsdebug !*arbvars !*varopt !*groebopt
        !*solveinconsistent depl!*);

%%% General changes

% Extend MAP/SELECT to other structures than list/matrix

symbolic procedure map!-eval1(o,q,fcn1,fcn2);
 % o       structure to be mapped.
 % q       map expression (univariate function).
 % fcn1    function for evaluating members of o.
 % fcn2    function computing results (e.g. aeval).
 map!-apply(map!-function(q,fcn1,fcn2),o);

symbolic procedure map!-function(q,fcn1,fcn2);
   begin scalar v,w;
   v := '!&!&x;
   if idp q 
      and (get(q,'simpfn) or get(q,'number!-of!-args)=1)
   then <<w:=v; q:={q,v}>> 
   else if eqcar(q,'replaceby) then
      <<w:=cadr q; q:=caddr q>>
   else
   <<w:=map!-frvarsof(q,nil);
      if null w then rederr "map/select: no free variable found" else
         if cdr w then rederr "map/select: free variable ambiguous";
      w := car w;
   >>;
   if eqcar(w,'!~) then w:=cadr w;
   q := sublis({w.v,{'!~,w}.v},q);
   return {'lambda,{'w},
      {'map!-eval2,'w,mkquote v,mkquote q,mkquote fcn1,mkquote fcn2}};
   end;

symbolic procedure map!-apply(f,o);
   if atom o then apply1(f,o)
   else (if m then apply2(m,f,o) 
      	 else car o . for each w in cdr o collect apply1(f,w))
       	    where m = get(car o,'mapfn);

symbolic procedure mapmat(f,o);
   'mat . for each row in cdr o collect
      for each w in row collect
	 apply1(f,w);

put('mat,'mapfn,'mapmat);


%%% EXCALC modifications

global '(indxl!*);
fluid '(kord!*);

% Remove covariant flag from indexed 0-forms.
% Add a subfunc to indexed forms.

if not getd 'excalcputform then copyd('excalcputform,'putform);

symbolic procedure putform(u,v);
   begin
   excalcputform(u,v);
   if not atom u then 
   << put(car u,'subfunc,'(lambda (a b) b));
      remflag({car u},'covariant) >>;
   end;


% Have to update ndepends to REDUCE3.6 depends.

symbolic procedure ndepends(u,v);
   if null u or numberp u or numberp v then nil
    else if u=v then u
    else if atom u and u memq frlis!* then t
      %to allow the most general pattern matching to occur;
    else if (lambda x; x and lndepends(cdr x,v)) assoc(u,depl!*)
     then t
    else if not atom u and idp car u and get(car u,'dname) then
        (if depends!-fn then apply2(depends!-fn,u,v) else nil)
           where (depends!-fn = get(car u,'domain!-depends!-fn))
    else if not atomf u
      and (lndepends(cdr u,v) or ndepends(car u,v)) then t
    else if atomf v or idp car v and get(car v,'dname) then nil
    % else ndependsl(u,cdr v);
    else nil;

%%% Make depends from ALG into ndepends from EXCALC (identical except
%%% for atomf test which treats indexed variables as atoms).

copyd('depends,'ndepends);

symbolic procedure lndepends(u,v);
   % Need to allow the possibility that U is an atom because the int
   % package calls depends with sq instead of prefix.
   if null u then nil
    else if atom u then ndepends(u,v)
    else ndepends(car u,v) or lndepends(cdr u,v);

% changed first u -> v (error otherwise)

symbolic procedure ndependsl(u,v);
   v and (ndepends(u,car v) or ndependsl(u,cdr v));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        OLD PATCHES: REMOVE ONCE IN PATCHES.RED!!!             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


fluid '(!*edsverbose !*edsdebug !*arbvars !*varopt !*groebopt
        !*solveinconsistent depl!*);


%%% SOLVE changes

% Changed depl!* line to work for non-atom kernels as well.

fluid '(!*expandexpt);   % from simp.red

fluid '( system!*        % system to be solved
         osystem!*       % original system on input
         uv!*            % user supplied variables
         iv!*            % internal variables
         fv!*            % restricted variables
         kl!*            % kernels to be investigated
         sub!*           % global substitutions
         inv!*           % global inverse substitutions
	 depl!*          % REDUCE dependency list
	 !*solvealgp     % true if using this module
         solvealgdb!*    % collecting some data
         last!-vars!*    % collection of innermost aux variables
         const!-vars!*   % variables representing constants
         root!-vars!*    % variables representing root expressions
         !*expli         % local switch: explicit solution
         groebroots!*    % predefined roots from input surds
         !*test_solvealg % debugging support
         !*arbvars
       );
 
fluid '(!*trnonlnr);
  % If set on, the modified system and the Groebner result
  % or the reason for the failure are printed. 

global '(loaded!-packages!* !!arbint);

symbolic procedure solvenonlnrsys2();
  % Main driver. We need non-local exits here
  % because of possibly hidden non algebraic variable
  % dependencies.
  if null !*solvealgp then system!*:='(FAILED) else % against recursion.
  (begin scalar iv!*,kl!*,inv!*,fv!*,r,w,!*solvealgp,solvealgdb!*,sub!*;
	 scalar last!-vars!*,groebroots!*,const!-vars!*,root!-vars!*;
         % preserving the variable sequence if *varopt is off
%      if not !*varopt then depl!* :=
%        append(pair(uv!*,append(cdr uv!*,{gensym()})),depl!*);
      if not !*varopt then depl!* :=
        append(foreach l on uv!* collect l,depl!*);
         % hiding dmode because exponentials need integers.
      for each f in system!* do solvealgk0 
         (if dmode!* then numr subf(f,nil) where dmode!*=nil else f);
      if !*trnonlnr then print list("original system:",system!*);
      if !*trnonlnr then print list("original kernels:",kl!*);
      if null cdr system!* then
          if (smemq('sin,system!*)or smemq('cos,system!*)) and
             (r:=solvenonlnrtansub(prepf(w:=car system!*),car uv!*))
             and car r
            then return solvenonlnrtansolve(r,car uv!*,w)
           else if (smemq('sinh,system!*)or smemq('cosh,system!*)) and
             (r:=solvenonlnrtanhsub(prepf(w:=car system!*),car uv!*))
             and car r
            then return solvenonlnrtanhsolve(r,car uv!*,w);
      if atom (errorset('(solvealgk1),!*trnonlnr,nil))
	       where dmode!*=nil
         then return (system!*:='(FAILED));
      system!*:='LIST.for each p in system!* collect prepf p;
      if not('groebner memq loaded!-packages!*)
        then load!-package 'groebner;
      for each x in iv!* do if not member(x,last!-vars!*) then
        for each y in last!-vars!* do depend1(x,y,t);
      iv!* := sort(iv!*,function (lambda(a,b);depends(a,b)));
      if !*trnonlnr then
      <<  prin2t "Entering Groebner for system";
          writepri(mkquote system!*,'only);
          writepri(mkquote('LIST.iv!*), 'only);
      >>;
      r := list(system!*,'LIST.iv!*);
      r := groesolveeval r;
      if !*trnonlnr then
      <<  prin2t "leaving Groebner with intermediate result";
          writepri(mkquote r,'only);
          terpri(); terpri();
      >>;
      if 'sin memq solvealgdb!* then r:=solvealgtrig2 r;
      if 'sinh memq solvealgdb!* then r:=solvealghyp2 r;
      r:= if r='(LIST) then '(INCONSISTENT) else solvealginv r;
      system!* := r;  % set value aside
      return r;
  end) where depl!*=depl!* ;


% Make variable dependency override reordering.


fluid '(!*trsparse);

symbolic procedure solvesparsecheck(sys,vl);
   % sys: list of sf, vl: list of kernel
   % -> solvesparsecheck: nil or {list of sf,list of kernel}
   % This program checks for a sparse linear system. If the
   % system is sparse enough, it returns (exlis.varlis) reordered
   % such that a maximum triangular upper diagonal form is
   % established. Otherwise the result is NIL.
   begin scalar vl1,xl,sys1,q,x,y;
         integer sp;    

   % First collect a list vl1 where each variable is followed
   % by the number of equations where it occurs, and then
   % by the number of other variables which occur in these
   % equations (connectivity). At the same time, collect a measure
   % of the sparsity.
   sp:=0;
   vl1:= for each x in vl collect x . 0 . nil;
   foreach q in sys do
      foreach x in (xl := intersection(topkerns q,vl)) do
       <<y := assoc(x,vl1);
         cdr y := (cadr y + 1) . union(xl,cddr y);
         sp := sp + 1>>;
   foreach p in vl1 do
      cddr p := length cddr p - 1; % could drop the -1

   % Drop out if density > 80%
   if sp > length sys * length vl * 0.8 then 
    <<if !*trsparse then prin2t "System is not very sparse";
      return nil>>;

   % If varopt is on, sort variables first by least occurrences and then
   % least connectivity, but allow dependency to override.
   % Reset kernel order and reorder equations.
   if !*trsparse then
     solvesparseprint("Original sparse system",sys,vl);
   
   if !*varopt then
   << vl1 := sort(vl1,function solvevarordp);
      vl1 := foreach x in vl1 collect car x;
%      if !*trsparse then lpriw("Optimal variable order:",vl1); 
%      foreach k in reverse vl1 do updkorder k;
%      vl1 := sort(vl1,function solvevarordp1);
      vl1 := solvevaradjust vl1;
%      if !*trsparse then lpriw("Adjusted variable order:",vl1); 
      foreach k in reverse vl1 do updkorder k;
      sys := for each q in sys collect reorder q >>
   else vl1 := foreach x in vl1 collect car x;

   % Next sort equations in ascending order of their first variable
   % and then descending order of the next variable.
   sys1:= (nil . nil) . foreach x in vl1 collect x . nil;
   foreach q in sys do
      <<if domainp q or not member(mvar q,vl1) then y := assoc(nil,sys1)
        else y := assoc(mvar q,sys1);
        cdr y := q . cdr y>>;
   foreach p in cdr sys1 do
      if cdr p then cdr p := sort(cdr p, function solvesparsesort);

   % Finally split off a leading diagonal system and push the remaining
   % equations down.
   sys := nconc(foreach p in sys1 join if cdr p then {cadr p},
                reversip foreach p in sys1 join if cdr p then cddr p);
   if !*trsparse then
      solvesparseprint("Variables and/or equations rearranged",sys,vl1);
   return sys.vl1;
   end;


symbolic procedure solvevarordp(x,y);
   cadr x < cadr y or
   cadr x = cadr y and cddr x < cddr y;

symbolic procedure solvevarordp1(x,y);
   % This is incomplete, since it is not transitive
   depends(x,y) or
   not depends(y,x) and ordop(x,y);


symbolic procedure solvevaradjust u;
   begin scalar v,y;
   if null u then return nil;
   v := foreach x in u join
      	<< y := assoc(x,depl!*);
 	   if null y or null xnp(cdr y,u) then {x} >>;
   return nconc(solvevaradjust setdiff(u,v),v);
   end;

% Usually solve goes to the Cramer method since the expressions
% contain exponentials. The Bareiss code should work, so disable this.

symbolic procedure exptexpflistp u;
   nil;

endmodule;

end;