reduce-historical/r38/packages/misc/pf.red

module pf;  % Compute partial fractions for an expression.

% Author: Anthony C. Hearn.

Comment PF is the top level operator for finding the partial fractions
of an expression.  It returns the partial fractions as a list.

The algorithms used here are relatively unsophisticated, and use the
extended Euclidean algorithm to break up expressions into factors.
Much more sophisticated algorithms exist in the literature;

fluid '(!*exp !*limitedfactors !*gcd kord!*);

symbolic operator pf;

flag('(pf),'noval);   % Since PF will do its own simplification.

symbolic procedure pf(u,var);
   % Convert an algebraic expression into partial fractions.
   begin scalar !*exp,!*gcd,kord!*,!*limitedfactors,polypart,rfactor,
		u1,u2,u3,u4,var,x,xx,y;
      !*exp := !*gcd := t;
      xx := updkorder var;         % Make var the main variable.
      x := subs2 resimp simp!* u;  % To allow for OFF EXP forms.
      u1 := denr x;
      if degr(u1,var) = 0 then <<setkorder xx; return list('list,u)>>;
      u2 := qremsq(!*f2q numr x,!*f2q u1,var); %Extract polynomial part.
      if caar u2 then polypart := car u2;
      rfactor := 1 ./ 1;           % Factor for rational part.
      u2 := cdr u2;
      u3 := fctrf u1;              % Factorize denominator.
      x := cdr u3;
      u3 := car u3;
      % Process monomial part.
      while not domainp u3 do
       <<if mvar u3 eq var then x := (!*k2f mvar u3 . ldeg u3) . x
	  else <<u4 := !*p2f lpow u3;
		 rfactor := numr rfactor ./ multf(u4,denr rfactor);
		 u1 := quotf(u1,u4)>>;
	 u3 := lc u3>>;
      if u3 neq 1 then <<rfactor := numr rfactor
				       ./ multf(u3,denr rfactor);
			 u1 := quotf(u1,u3)>>;
      % Separate power factors in denominator.
      while length x>1 do
       <<u3 := exptf(caar x,cdar x);
	 u1 := quotf(u1,u3);
	 if degr(u3,var)=0
	   then rfactor := numr rfactor ./ multf(u3,denr rfactor)
      %    then <<rfactor := numr rfactor ./ multf(u3,denr rfactor);
      %           u2 := nil>>
	  else <<u4 := xeucl(u1,u3,var);
		 % Remove spurious polynomial in numerator.
		 y := (multsq(remsq(multsq(car u4,u2),!*f2q u3,var),
			      rfactor) . car x)
			. y;
		 u2 := multsq(cdr u4,u2)>>;
	 x := cdr x>>;
      u3 := exptf(caar x,cdar x);
      if u2 = (nil ./ 1) then nil
       else if degr(u3,var)=0
	then rfactor := numr rfactor ./ multf(u3,denr rfactor)
      % Remove spurious polynomial in numerator.
       else y := (multsq(rfactor,remsq(u2,!*f2q u3,var)) . car x) . y;
      x := nil;
      % Finally break down non-linear terms in denominator.
      for each j in y do
	 if cddr j =1 then x := j . x
	  else x := append(pfpower(car j,cadr j,cddr j,var),x);
      x := for each j in x
	      collect list('quotient,prepsq!* car j,
			   if cddr j=1 then prepf cadr j
			    else list('expt,prepf cadr j,cddr j));
      if polypart then x := prepsq!* polypart . x;
      setkorder xx;
      return 'list . x
   end;

symbolic procedure xeucl(u,v,var);
   % Extended Euclidean algorithm with rational coefficients.
   % I.e., find polynomials Q, R in var with rational coefficients (as 
   % standard quotients) such that Q*u + R*v = 1, where u and v are
   % relatively prime standard forms in variable var.  Returns Q . R.
   begin scalar q,r,s,w;
      q := list(1 ./ 1,nil ./ 1);
      r := list(nil ./ 1,1 ./ 1);
      if degr(u,var) < degr(v,var)
	then <<s := u; u := v; v := s; s := q; q := r; r := s>>;
      u := !*f2q u; v := !*f2q v;
      while numr v do
       <<if degr(numr v,var)=0 then w := quotsq(u,v) . (nil ./ 1)
	  else w := qremsq(u,v,var);
	 s := list(addsq(car q,negsq multsq(car w,car r)),
		   addsq(cadr q,negsq multsq(car w,cadr r)));
	 u := v;
	 v := cdr w;
	 q := r;
	 r := s>>;
      v := lnc numr u ./ denr u;   % Is it possible for this not to be
				   % in lowest terms, and, if so, does
				   % it matter?
      r := quotsq(v,u);
      return multsq(r,quotsq(car q,v)) . multsq(r,quotsq(cadr q,v))
  end;

symbolic procedure qremsq(u,v,var);
   % Find rational quotient and remainder (as standard quotients)
   % dividing standard quotients u by v wrt var.
   % This should really be done more directly without using quotsq.
   (quotsq(addsq(u,negsq x),v) . x) where x=remsq(u,v,var);

symbolic procedure remsq(u,v,var);
   % Find rational and remainder (as a standard quotient) on
   % dividing standard quotients u by v wrt var.
   begin integer m,n; scalar x;
      n := degr(numr v,var);
      if n=0 then rederr list "Remsq given zero degree polynomial";
      while (m := degr(numr u,var))>= n do
       <<if m=n then x := v
	  else x := multsq(!*p2q(var.**(m-n)),v);
	 u := addsq(u,
		    negsq multsq(multf(lc numr u,denr v)
				   ./ multf(lc numr v,denr u),
				 x))>>;
      return u
   end;

symbolic procedure pfpower(u,v,n,var);
   % Convert u/v^n into partial fractions.
   begin scalar x,z;
      while degr(numr u,var)>0 do
       <<x := qremsq(u,!*f2q v,var);
	 z := (cdr x . v . n) . z;
	 n := n-1;
	 u := car x>>;
      if numr u then z := (u . v . n) . z;
      return z
   end;

endmodule;

end;