module divide; % Exact division of standard forms to give a S Q.
% Authors: Mary Ann Moore and Arthur C. Norman.
fluid '(!*trdiv !*trint resid sqrtlist zlist);
exports fquotf,testdivdf,dfquotdf;
imports df2q,f2df,gcdf,interr,multdf,negdf,plusdf,printdf,printsf,
quotf,multsq,invsq,negsq;
% Intended for dividing out known factors as produced by the
% integration program. Horrible and slow, I expect!!
symbolic procedure dfquotdf(a,b);
begin scalar resid;
if (!*trint or !*trdiv) then <<
printc "Dividing out a simple factor of "; printdf b >>;
a:=dfquotdf1(a,b);
if (!*trint or !*trdiv) then <<
printc "Remaining term to be factorised is ";
printdf a >>;
if not null resid then begin
scalar gres,w;
% Make one more check for a zero residue.
if null numr df2q resid then return nil;
if !*trint or !*trdiv then <<
printc "Failure in factorisation:";
printdf resid;
printc "Which should be zero";
w:=resid;
gres:=numr lc w; w:=red w;
while not null w do <<
gres:=gcdf(gres,numr lc w);
w:=red w >>;
printc "I.e. the following vanishes";
printsf gres>>;
interr "Non-exact division due to a log term"
end;
return a
end;
symbolic procedure fquotf(a,b);
% Input: a and b standard quotients with (a/b) an exact
% division with respect to the variables in zlist,
% but not necessarily obviously so. the 'non-obvious' problems
% will be because of (e.g.) square-root symbols in b
% output: standard quotient for (a/b)
% (prints message if remainder is not 'clearly' zero.
% A must not be zero.
begin scalar t1;
if null a then interr "a=0 in fquotf";
t1:=quotf(a,b); %try it the easy way
if not null t1 then return t1 ./ 1; %ok
return df2q dfquotdf(f2df a,f2df b)
end;
symbolic procedure dfquotdf1(a,b);
begin scalar q;
if null b then interr "Attempt to divide by zero";
q:=sqrtlist; %remove sqrts from denominator, maybe.
while not null q do begin
scalar conj;
conj:=conjsqrt(b,car q); %conjugate wrt given sqrt
if not (b=conj) then <<
a:=multdf(a,conj);
b:=multdf(b,conj) >>;
q:=cdr q end;
q:=dfquotdf2(a,b);
resid:=reversip resid;
return q
end;
symbolic procedure dfquotdf2(a,b);
% As above but a and b are distributed forms, as is the result.
if null a then nil
else begin scalar xd,lcd;
xd:=xpdiff(lpow a,lpow b);
if xd='failed then <<
xd:=lt a; a:=red a;
resid:=xd .+ resid;
return dfquotdf2(a,b) >>;
lcd:= !*multsq(lc a,!*invsq lc b);
if null numr lcd then return dfquotdf2(red a,b);
% Should not be necessary;
lcd := xd .* lcd;
xd:=plusdf(a,multdf(negdf (lcd .+ nil),b));
if xd and (lpow xd = lpow a % Again, should not be necessary;
or xpdiff(lpow xd,lpow b) = 'failed)
then <<if !*trint or !*trdiv
then <<printc "Problems in dividing:"; printdf xd>>;
xd := rootextractdf xd;
if !*trint or !*trdiv then printdf xd>>;
return lcd .+ dfquotdf2(xd,b)
end;
symbolic procedure rootextractdf u;
if null u then nil
else begin scalar v;
v := resimp rootextractsq lc u;
return if null numr v then rootextractdf red u
else (lpow u .* v) .+ rootextractdf red u
end;
symbolic procedure rootextractsq u;
if null numr u then u
% else rootextractf numr u ./ rootextractf denr u;
else (rootextractf numr x ./ rootextractf denr x)
where x=subs2q u;
symbolic procedure rootextractf v;
if domainp v then v
else begin scalar u,r,c,x,p;
u := mvar v; p := ldeg v;
r := rootextractf red v;
c := rootextractf lc v;
if null c then return r
else if atom u then return (lpow v .* c) .+ r
else if car u eq 'sqrt
or car u eq 'expt and eqcar(caddr u,'quotient)
and car cdaddr u = 1 and numberp cadr cdaddr u
then <<p := divide(p,if car u eq 'sqrt then 2
else cadr cdaddr u);
if car p = 0
then return if null c then r else (lpow v .* c) .+ r
else if numberp cadr u
then <<c := multd(cadr u ** car p,c); p := cdr p>>
else <<x := simpexpt list(cadr u,car p);
if denr x = 1 then <<c := multf(numr x,c); p := cdr p>>
else p := ldeg v>>>>;
% D. Dahm suggested an additional call of rootextractf on the
% result here. This does cause some expressions to simplify
% sooner, but also leads to infinite loops with expressions
% like (a*x+b)**p.
return if p=0 then addf(c,r)
else if null c then r
else ((u to p) .* c) .+ r
end;
% The following hack makes sure that the results of differentiation
% gets passed through ROOTEXTRACT
% a) This should not be done this way, since the effect is global
% b) Should this be done via TIDYSQRT?
put('df,'simpfn,'simpdf!*);
symbolic procedure simpdf!* u;
begin scalar v,v1;
v:=simpdf u;
v1:=rootextractsq v;
if not(v1=v) then return resimp v1
else return v
end;
symbolic procedure xpdiff(a,b);
%Result is list a-b, or 'failed' if a member of this would be negative.
if null a then if null b then nil
else interr "B too long in xpdiff"
else if null b then interr "A too long in xpdiff"
else if car b>car a then 'failed
else (lambda r;
if r='failed then 'failed
else (car a-car b) . r) (xpdiff(cdr a,cdr b));
symbolic procedure conjsqrt(b,var);
% Subst(var=-var,b).
if null b then nil
else conjterm(lpow b,lc b,var) .+ conjsqrt(red b,var);
symbolic procedure conjterm(xl,coef,var);
% Ditto but working on a term.
if involvesp(xl,var,zlist) then xl .* negsq coef
else xl .* coef;
symbolic procedure involvesp(xl,var,zl);
% Check if exponent list has non-zero power for variable.
if null xl then interr "Var not found in involvesp"
else if car zl=var then car xl neq 0
else involvesp(cdr xl,var,cdr zl);
endmodule;
end;