module intfac; % Interface between integrator and factorizer.
% Author: Anthony C. Hearn.
% Based on earlier versions by James Davenport, Mary Ann Moore and
% Arthur Norman.
fluid '(!*intfac !*surds kord!* zlist); % clogflag
exports int!-fac;
symbolic procedure int!-fac x;
% X is a primitive, square-free polynomial, except for monomial
% factors. Result is a list of 'factors' wrt zlist.
% Throughout most of the integrator we want to add new surds, so we
% turn surds on. However, we use *intfac to inhibit use of quadratic
% factorizer in the poly/primfac module, since things don't work
% properly if this is used.
begin scalar !*intfac,!*surds;
!*intfac := !*surds := t;
return int!-fac!-inner x;
end;
symbolic procedure int!-fac!-inner x;
% X is a primitive, square-free polynomial, except for monomial
% factors. Result is a list of 'factors' wrt zlist.
begin scalar factors;
factors := fctrf x;
factors := cdr factors; % Ignore monomial factor.
% Make sure x really square-free.
factors := for each u in factors
collect if cdr u=1 then car u
else interr list(x,"not square free");
% It seems we need the logs ordered ahead of atans, hence reverse.
return reversip for each u in factors join fac2int u
end;
symbolic procedure fac2int u;
% Returns a list of all the arctangents and logarithms arising from
% an attempt to take the one irreducible (but not necessarily the
% absolutely irreducible) factor u.
begin scalar degrees,x;
degrees := for each w in zlist collect (degreef(u,w) . w);
if assoc(1,degrees) then return list ('log . (u ./ 1))
% An irreducible polynomial of degree 1 is absolutely irreducible.
else if x := assoc(2,degrees) then return int!-quadterm(u,cdr x)
else if assoc(0,degrees) then return list('log . (u ./ 1));
% This suggests a surd occurs. Should that be an error?
if !*trint
then <<printc "*** Polynomial";
printsf u;
printc "has not been completely factored">>;
return list ('log . (u ./ 1))
end;
symbolic procedure int!-quadterm(pol,var);
% Add in logs and atans corresponding to splitting the polynomial pol
% given it is quadratic wrt var. Does not assume pol is univariate.
% We need to rootxf!* so that
% int(1/(x**2*y0+x**2+x*y0**2+3*x*y0+x+y0**2+y0),x) comes out in
% terms of real functions.
begin scalar a,b,c,discrim,kord,res,w;
kord := setkorder(var . kord!*); % It shouldn't matter if
% var occurs twice.
c := reorder pol;
if ldeg c neq 2
then <<setkorder kord;
rerror(int,5,"Invalid polynomial in int-quadterm")>>;
a := lc c;
c := red c;
if not domainp c and mvar c = var and ldeg c = 1
then <<b := lc c; c := red c>>;
setkorder kord;
discrim := powsubsf addf(multf(b,b),multd(-4,multf(a,c)));
if null discrim then interr "discrim is zero in quadterm";
% A quadratic usually implies an atan term.
% if not clogflag
% then <<w := rootxf(negf discrim,2);
% if not(w eq 'failed) then go to atancase>>;
w := rootxf!*(negf discrim,2);
if not(w eq 'failed) then go to atancase;
w := rootxf!*(discrim,2); % Maybe only rootxf is needed here.
if w eq 'failed then return list ('log . !*f2q pol);
% if w eq 'failed then rederr "Integration failure in int-quadterm";
discrim := w;
w := multpf(mksp(var,1),a);
w := addf(multd(2,w),b); % 2*a*x + b.
a := addf(w,discrim);
b := addf(w,negf discrim);
% Remove monomial multipliers.
a := quotf(a,cdr comfac a);
b := quotf(b,cdr comfac b);
return ('log . !*f2q a) . ('log . !*f2q b) . res;
atancase:
res := ('log . !*f2q pol) . res; % One part of answer.
a := multpf(mksp(var,1),a);
a := addf(b,multd(2,a));
a := fquotf(a,w);
return ('atan . a) . res
end;
symbolic procedure rootxf!*(u,n);
(if x eq 'failed or smemq('i,x) and not smemq('i,u)
then (rootxf(u,n) where !*surds=nil)
else x)
where x=rootxf(u,n);
endmodule;
end;