reduce-historical/r37/packages/matrix/resultnt.red

module resultnt;

% Author: Eberhard Schruefer.

% Modifications by: Anthony C. Hearn, Winfried Neun.

%**********************************************************************
%                                                                     *
% The resultant function defined here has the following properties:   *
%                                                                     *
%                           degr(p1,x)*degr(p2,x)                     *
%  resultant(p1,p2,x) = (-1)                     *resultant(p2,p1,x)  *
%                                                                     *
%                         degr(p2,x)                                  *
%  resultant(p1,p2,x) = p1             if p1 free of x                *
%                                                                     *
%  resultant(p1,p2,x) = 1  if p1 free of x and p2 free of x           *
%                                                                     *
%**********************************************************************

%exports resultant;

%imports reorder,setkorder,degr,addf,negf,multf,multpf;

load_package polydiv;

fluid '(!*bezout !*exp kord!*);

switch bezout;

put('resultant,'simpfn,'simpresultant);

symbolic procedure simpresultant u;
   if length u neq 3
     then rerror(matrix,19,
		   "Resultant called with wrong number of arguments")
    else resultantsq(simp!* car u,simp!* cadr u,!*a2k caddr u)
	 where !*exp = t;

symbolic procedure resultant(u,v,var);
   % Kept for compatibility with old code.
   if domainp u and domainp v then 1
    else begin scalar x;
       kord!* := var . kord!*;  % updkorder can't be used here.
				% See sum test.
       if null domainp u and null(mvar u eq var) then u := reorder u;
       if null domainp v and null(mvar v eq var) then v := reorder v;
       x := if !*bezout then bezout_resultant(u,v,var)
	     else !*a2f polyresultant(prepf u,prepf v,var);
       setkorder cdr kord!*;
       return x
   end;

symbolic procedure resultantsq(u,v,var);
   if domainp numr u and domainp numr v and denr u = 1 and denr v = 1
     then 1 ./ 1
    else begin scalar x;
       kord!* := var . kord!*;  % updkorder can't be used here.
				% See sum test.
       if null domainp numr u and null(mvar numr u eq var)
	 then u := reordsq u;
       if null domainp numr v and null(mvar numr v eq var)
	 then v := reordsq v;
       x := if !*bezout then !*f2q bezout_resultant(!*q2f u,!*q2f v,var)
	     else simp polyresultant(prepsq u,prepsq v,var);
       setkorder cdr kord!*;
       return x
   end;

algebraic (co_off := { co(0,~x) => x });

% algebraic procedure notunivariatep(uu);
%    for i:=1:length uu sum if fixp part(uu,i) then 0 else 1;
   
algebraic procedure notunivariatep uu;
   for each u in uu sum if fixp u then 0 else 1;

algebraic procedure polyresultant(u,v,var);
   % See Zippel's book p 151, subresultant algorithm --
   % more or less the same.
   begin scalar g,h,delta,beta,temp,uu,vv;
      uu := coeff(u,var); vv := coeff(v,var);
      if length uu<length vv
	then return (-1 * polyresultant(v,u,var))
       else if (notunivariatep(uu) > 0) or (notunivariatep(vv)>0)
	then <<u := for i:=1:length uu sum
			(if fixp part(uu,i) then part(uu,i)
			  else (co(0,part(uu,i))))*var^(i-1);
	       v := for i:=1:length vv sum
			(if fixp part(vv,i) then part(vv,i)
			  else (co(0,part(vv,i))))*var^(i-1)>>;
      % Change to nested domain.
      g := h := 1;
      while not (V=0) do
       <<delta := deg(u,var) - deg(v,var);
	 beta := (-1)^(delta +1) * g * h^delta;
	 h := h*(lcof(v,var)/h)^delta;
	 temp := u;
	 u := v;
	 v := pseudo_remainder(temp,v,var)/beta;
	 g := lcof(u,var)>>;
	 if deg(u,var) = 0 then u := u^delta else return 0;
      let co_off; u := u; clearrules co_off;
      return u
   end;

symbolic procedure lcoff(u,var);
   if domainp u or not(mvar u eq var) then 0 else lc u;

symbolic procedure bezout_resultant(u,v,w);
   % U and v are standard forms. Result is resultant of u and v
   % w.r.t. kernel w. Method is Bezout's determinant using exterior
   % multiplication for its calculation.
   begin integer n,nm; scalar ap,ep,uh,ut,vh,vt;
      if domainp u or null(mvar u eq w)
	then return if not domainp v and mvar v eq w 
                      then exptf(u,ldeg v)
		     else 1
       else if domainp v or null(mvar v eq w)
	then return if mvar u eq w then exptf(v,ldeg u) else 1;
      n := ldeg u - ldeg v;
      ep := 1;
      if n<0 then
          <<for j := (-n-1) step -1 until 1 do
              ep := b!:extmult(!*sf2exb(multpf(w to j,u),w),ep);
	      ep := b!:extmult(!*sf2exb(multd((-1)**(-n*ldeg u),u),
					w),
			       ep)>>
       else if n>0 then
            <<for j := (n-1) step -1 until 1 do
                ep := b!:extmult(!*sf2exb(multpf(w to j,v),w),ep);
              ep := b!:extmult(!*sf2exb(v,w),ep)>>;
     nm := max(ldeg u,ldeg v);
     uh := lc u;
     vh := lc v;
     ut := if n<0 then multpf(w to -n,red u)
	   else red u;
     vt := if n>0 then multpf(w to n,red v)
            else red v;
     ap := addf(multf(uh,vt),negf multf(vh,ut));
     ep := if null ep then !*sf2exb(ap,w)
        else b!:extmult(!*sf2exb(ap,w),ep);
     for j := (nm - 1) step -1 until (abs n + 1) do
        <<if degr(ut,w) = j then
         <<uh := addf(lc ut,multf(!*k2f w,uh));
                   ut := red ut>>
       else    uh := multf(!*k2f w,uh);
          if degr(vt,w) = j then
         <<vh := addf(lc vt,multf(!*k2f w,vh));
                   vt := red vt>>
       else    vh := multf(!*k2f w,vh);
      ep := b!:extmult(!*sf2exb(addf(multf(uh,vt),
                    negf multf(vh,ut)),w),ep)>>;
      return if null ep then nil else lc ep
   end;

symbolic procedure !*sf2exb(u,v);
   %distributes s.f. u with respect to powers in v.
   if degr(u,v)=0 then if null u then nil
                        else list 0 .* u .+ nil
    else list ldeg u .* lc u .+ !*sf2exb(red u,v);

%**** Support for exterior multiplication ****
% Data structure is lpow ::= list of degrees in exterior product
%                   lc   ::= standard form

symbolic procedure b!:extmult(u,v);
   %Special exterior multiplication routine. Degree of form v is
   %arbitrary, u is a one-form.
   if null u or null v then  nil
    else if v = 1 then u
    else (if x then cdr x .* (if car x then negf multf(lc u,lc v)
                   else multf(lc u,lc v))
              .+ b!:extadd(b!:extmult(!*t2f lt u,red v),
                    b!:extmult(red u,v))
       else b!:extadd(b!:extmult(red u,v),
              b!:extmult(!*t2f lt u,red v)))
      where x = b!:ordexn(car lpow u,lpow v);

symbolic procedure b!:extadd(u,v);
   if null u then v
    else if null v then u
    else if lpow u = lpow v then
            (lambda x,y; if null x then y else lpow u .* x .+ y)
        (addf(lc u,lc v),b!:extadd(red u,red v))
    else if b!:ordexp(lpow u,lpow v) then lt u .+ b!:extadd(red u,v)
    else lt v .+ b!:extadd(u,red v);

symbolic procedure b!:ordexp(u,v);
   if null u then t
    else if car u > car v then t
    else if car u = car v then b!:ordexp(cdr u,cdr v)
    else nil;

symbolic procedure b!:ordexn(u,v);
   %u is a single integer, v a list. Returns nil if u is a member
   %of v or a dotted pair of a permutation indicator and the ordered
   %list of u merged into v.
   begin scalar s,x;
     a: if null v then return(s . reverse(u . x))
     else if u = car v then return nil
     else if u and u > car v then
                 return(s . append(reverse(u . x),v))
         else  <<x := car v . x;
                 v := cdr v;
                 s := not s>>;
         go to a
   end;

endmodule;

end;