module solve;   % Solve one or more algebraic equations.

% Author: David R. Stoutemyer.
% Major modifications by: David Hartley, Anthony C. Hearn, Herbert
% Melenk, Donald R. Morrison and Rainer Schoepf.

create!-package('(solve solve1 ppsoln solvelnr glsolve solvealg solvetab
		  quartic),nil);

% Other packages needed by solve package.

load!-package 'matrix;

fluid '(!*exp !*ezgcd !*multiplicities !!gcd dmode!* vars!*);

fluid '(inside!-solveeval solve!-gensymcounter);

% solve!-gensymcounter := 1;

global '(multiplicities!* assumptions requirements);

flag('(multiplicities!* assumptions requirements),
     'share);

% Those switches that are on are now set in entry.red.

% !*multiplicities      Lists all roots with multiplicities if on.
%  !!gcd                SOLVECOEFF returns GCD of powers of its arg in
%                       this.  With the decompose code, this should
%                       only occur with expressions of form x^n + c.

algebraic operator one_of;

put('arbint,'simpfn,'simpiden);

% algebraic operator arbreal;

symbolic operator expand_cases;

symbolic procedure simp!-arbcomplex u;
    simpiden('arbcomplex . u) where dmode!*=nil;

deflist('((arbcomplex simp!-arbcomplex)),'simpfn);


% ***** Utility Functions *****

symbolic procedure freeofl(u,v);
   null v or freeof(u,car v) and freeofl(u,cdr v);

symbolic procedure allkern elst;
   % Returns list of all top-level kernels in the list of standard
   % quotients elst.   Corrected 5 Feb 92 by Francis Wright.
   if null elst then nil
    else union(kernels numr car elst, allkern cdr elst);

symbolic procedure topkern(u,x);
   % Returns list of top level kernels in the standard form u that
   % contain the kernel x;
   for each j in kernels u conc if not freeof(j,x) then list j else nil;

symbolic procedure coeflis ex;
   % Ex is a standard form.  Returns a list of the coefficients of the
   % main variable in ex in the form ((expon . coeff) (expon . coeff)
   % ... ), where the expon's occur in increasing order, and entries do
   % not occur of zero coefficients.  We need to reorder coefficients
   % since kernel order can change in the calling function.
   begin scalar ans,var;
      if domainp ex then return (0 . ex);
      var := mvar ex;
      while not domainp ex and mvar ex=var do
	<<ans := (ldeg ex . reorder lc ex) . ans; ex := red ex>>;
      if ex then ans := (0 . reorder ex) . ans;
      return ans
   end;


% ***** Evaluation Interface *****

% Solvemethods!* is a list of procedures which are able to process
% one problem class. Each of its members must check itself
% whether it can be applied or not. The classical equation solver
% is called if none of the methods can contribute.
%
% Protocol:
%
%   input: PSOPFN standard, where the elements of the input list
%          have been passed through REVAL.
%
%   output:
%          'nil: the algorithm cannot be applied because the problem
%               belongs to a different problem class;
%          '(failed): the problem belongs to the class represented
%               by the algorithm but the program has been
%               unable to compute a result. The problem should
%               not be given to any other method - instead the
%               input should be returned.
%          result: the algorithm has been successful and the final
%               result is returned as algebraic form (including an
%               eventually empty result for an "inconsistent" case).

fluid '(solvemethods!*);

put('solve,'psopfn,'solveeval);

symbolic procedure solveeval u;
  begin scalar w,r,m;
    w:=for each q in u collect reval q;
    m:=solvemethods!*;
    while null r and m do <<r := apply1(car m,w); m := cdr m>>;
    return if null r then solveeval1 w
      else if eqcar(r,'failed) then 'solve . u
      else r;
  end;

% Links to other packages.

symbolic procedure odesolve!* u;
   % 2 arg solve => algebraic always (otherwise cannot algebraically
   % solve an equation involving a derivative!)
   length u neq 2 and smemq('df,u) and
   <<load!-package 'odesolve;
     % Support both "old" ODESOLVE and ODESolve 1.04+:
     if flagp('odesolve, 'opfn) and length u neq 3 then '(failed)
      else aeval('odesolve . u)>>;

solvemethods!* := union('(odesolve!*),solvemethods!*);

symbolic procedure solveeval1 u;
   begin scalar !*ezgcd,!!gcd,vars!*;  integer nargs;
      if atom u then rerror(solve,1,"SOLVE called with no equations")
       else if null dmode!* then !*ezgcd := t;
      nargs := length u;
      if not inside!-solveeval then 
      <<solve!-gensymcounter := 1; 
        assumptions :=requirements:={'list}>>;
      u := (if nargs=1 then solve0(car u,nil)
              else if nargs=2 then solve0(car u, cadr u)
              else <<lprim "Please put SOLVE unknowns in a list";
                     solve0(car u,'list . cdr u)>>)
             where inside!-solveeval = t, !*resimp = not !*exp;
      if not inside!-solveeval then
      <<assumptions := solve!-clean!-info(assumptions,t);
        requirements:= solve!-clean!-info(requirements,nil)>>;
      return !*solvelist2solveeqlist u
    end;

symbolic procedure !*solvelist2solveeqlist u;
   begin scalar x,y,z;
      u := for each j in u collect solveorder j;
      for each j in u do
         <<if caddr j=0 then rerror(solve,2,"zero multiplicity")
            else if null cadr j
             then  x := for each k in car j collect
                                               list('equal,!*q2a k,0)
            else x := for each k in pair(cadr j,car j)
                          collect list('equal,car k,!*q2a cdr k);
           if length vars!* > 1 then x := 'list . x else x := car x;
           z := (caddr j . x) . z>>;
      z := sort(z,function ordp);
      x := nil;
      if !*multiplicities
	 then <<for each k in z do for i := 1:car k do x := cdr k . x;
		multiplicities!* := nil>>
       else <<for each k in z do << x := cdr k . x; y := car k . y>>;
	      multiplicities!* := 'list . reversip y>>;
      % Now check for redundant solutions.
%     if length vars!*>1 then z := check_solve_redundancy z;
      return 'list . reversip x
   end;

symbolic procedure solveorder u;
   % Put solve solutions in same order as specified variables.
   begin scalar v,w,x,y,z;
      v := vars!*;
      x := cadr u;      % SOLVE variable order.
      % Check if there are less variables than specified.
      if length x<length v then v := setdiff(v,setdiff(v,x));
      if null x or x = v then return u;
      y := car u;       % List of solutions.
      while x do <<z := (car x . car y) . z; x := cdr x; y := cdr y>>;
      w := v;
  a:  if null w then return reversip x . v . cddr u
       else if null(y := depassoc(car w,z)) then return u
       else x := cdr y . x;
      w := cdr w;
      go to a
   end;

symbolic procedure depassoc(u,v);
   if null v then nil
    else if u = caar v then car v
    else if depends(caar v,u) then nil   % Can't change order.
    else depassoc(u,cdr v);

% symbolic procedure check_solve_redundancy u;
%     % We assume all solutions are prefixed by LIST.
%     begin scalar x,y;
%        x := for each j in u collect cdr j;   %  Remove the LIST.
%        for each j in u do if not supersetlist(cdr j,x) then y:= j . y;
%        return reversip!* y
%     end;

% symbolic procedure supersetlist(u,v);
%    % Returns true if u is a non-equal superset of any element of v.
%    v and
%      (u neq car v and null setdiff(car v,u) or supersetlist(u,cdr v));

endmodule;

end;
