module int!-intro; % General support for REDUCE integrator.
% Authors: A. C. Norman and P. M. A. Moore.
% Modified by: J. Davenport, J. P. Fitch, A. C. Hearn.
% Note that at one point, INT had been flagged SIMP0FN. However, that
% lead to problems when the arguments of INT contained pattern
% variables.
fluid '(!*conscount !*noextend !*pvar gaussiani);
global '(btrlevel frlis!* gensymcount initl!*);
!*conscount:=10000; % default maximum number of conses in certain
% operations.
!*pvar:='!_a;
btrlevel := 5; %default to a reasonably full backtrace.
% The next smacro is needed at this point to define gaussiani.
symbolic smacro procedure !*kk2f u; !*p2f mksp(u,1);
gaussiani := !*kk2f '(sqrt -1);
gensymcount := 0;
initl!* := append('(!*noextend), initl!*);
flag('(interr),'transfer); %For the compiler;
flag ('(atan dilog erf expint expt log tan),'transcendental);
comment Kludge to define derivative of an integral and integral of
a derivative;
frlis!* := union('(!=x !=y),frlis!*);
put('df,'opmtch,'(((int !=y !=x) !=x) (nil . t)
(evl!* !=y) nil) . get('df,'opmtch));
put('int,'opmtch,'(((df !=y !=x) !=x) (nil . t)
(evl!* !=y) nil) . get('int,'opmtch));
put('evl!*,'opmtch,'(((!=x) (nil . t) !=x nil)));
put('evl!*,'simpfn,'simpiden);
%Various functions used throughout the integrator.
smacro procedure !*kk2q a; ((mksp(a,1) .* 1) .+ nil) ./ 1;
symbolic smacro procedure divsf(u,v); sqrt2top(u ./ v);
symbolic procedure flatten u;
if null u then nil
else if atom u then list u
else if atom car u then car u . flatten cdr u
else nconc(flatten car u,flatten cdr u);
symbolic procedure int!-gensym1 u;
<< gensymcount:=gensymcount+1;
compress append(explode u,explode gensymcount) >>;
symbolic smacro procedure maninp(u,v,w);
interr "MANINP called -- not implemented";
symbolic procedure mknill n; if n=0 then nil else nil . mknill(n-1);
% Various selectors written as macros.
smacro procedure argof u;
% Argument of a unary function.
cadr u;
smacro procedure firstsubs u;
% The first substitution in a substitution list.
car u;
smacro procedure lsubs u; car u;
smacro procedure rsubs u; cdr u;
smacro procedure lfirstsubs u; caar u;
smacro procedure rfirstsubs u; cdar u;
put('nthroot,'simpfn,'simpiden);
% The binary n-th root operator nthroot(x,2)=sqrt(x)
% no simplification is used here.
% Hope is that pbuild introduces it, and simplog removes it.
% Selectors for the taylor series structure.
% Format is:
%function.((first.last computed so far) . assoc list of computed terms).
% ***store-hack-1***:
% remove this macro if more store is available.
smacro procedure tayshorten u;nil;
smacro procedure taylordefn u; car u;
symbolic smacro procedure taylorfunction u; caar u;
smacro procedure taylornumbers u; cadr u;
smacro procedure taylorfirst u; caadr u;
smacro procedure taylorlast u; cdadr u;
smacro procedure taylorlist u; cddr u;
smacro procedure taylormake(fn,nums,alist); fn.(nums.alist);
endmodule;
module contents;
% Authors: Mary Ann Moore and Arthur C. Norman
fluid '(clogflag content indexlist sqfr varlist zlist);
exports contents,contentsmv,dfnumr,difflogs,factorlistlist,multsqfree,
multup,sqfree,sqmerge;
imports int!-fac,fquotf,gcdf,interr,!*multf,partialdiff,quotf,ordop,
addf,negf,domainp,difff,mksp,negsq,invsq,addsq,!*multsq,diffsq;
comment we assume no power substitution is necessary in this module;
symbolic procedure contents(p,v);
% Find the contents of the polynomial p wrt variable v;
% Note that v may not be the main variable of p;
if domainp(p) then p
else if v=mvar p then contentsmv(p,v,nil)
else if ordop(v,mvar p) then p
else contentsmv(makemainvar(p,v),v,nil);
symbolic procedure contentsmv(p,v,sofar);
% Find contents of polynomial P;
% V is main variable of P;
% SOFAR is partial result;
if sofar=1 then 1
else if domainp p then gcdf(p,sofar)
else if not v=mvar p then gcdf(p,sofar)
else contentsmv(red p,v,gcdf(lc p,sofar));
symbolic procedure makemainvar(p,v);
% Bring v up to be the main variable in polynomial p;
% Note that the reconstructed p must be used with care since;
% it does not conform to the normal reduce ordering rules;
if domainp p then p
else if v=mvar p then p
else mergeadd(mulcoeffsby(makemainvar(lc p,v),lpow p,v),
makemainvar(red p,v),v);
symbolic procedure mulcoeffsby(p,pow,v);
% Multiply each coefficient in p by the standard power pow;
if null p then nil
else if domainp p or not v=mvar p then ((pow .* p) .+ nil)
else (lpow p .* ((pow .* lc p) .+ nil)) .+ mulcoeffsby(red p,pow,v);
symbolic procedure mergeadd(a,b,v);
% Add polynomials a and b given that they have same main variable v;
if domainp a or not v=mvar a then
if domainp b or not v=mvar b then addf(a,b)
else lt b .+ mergeadd(a,red b,v)
else if domainp b or not v=mvar b then
lt a .+ mergeadd(red a,b,v)
else (lambda xc;
if xc=0 then (lpow a .* addf(lc a,lc b)) .+
mergeadd(red a,red b,v)
else if xc>0 then lt a .+ mergeadd(red a,b,v)
else lt b .+ mergeadd(a,red b,v))
(tdeg lt a-tdeg lt b);
symbolic procedure sqfree(p,vl);
if (null vl) or (domainp p) then
<<content:=p; nil>>
else begin scalar w,v,dp,gg,pg,dpg,p1,w1;
w:=contents(p,car vl); % content of p ;
p:=quotf(p,w); % make p primitive;
w:=sqfree(w,cdr vl); % process content by recursion;
if p=1 then return w;
v:=car vl; % pick out variable from list;
while not (p=1) do <<
dp:=partialdiff(p,v);
gg:=gcdf(p,dp);
pg:=quotf(p,gg);
dpg:=negf partialdiff(pg,v);
p1:=gcdf(pg,addf(quotf(dp,gg),dpg));
w1:=p1.w1;
p:=gg>>;
return sqmerge(reverse w1,w,t)
end;
symbolic procedure sqmerge(w1,w,simplew1);
% w and w1 are lists of factors of each power. if simplew1 is true
% then w1 contains only single factors for each power. ;
if null w1 then w
else if null w then if car w1=1 then nil.sqmerge(cdr w1,w,simplew1)
else (if simplew1 then list car w1 else car w1).
sqmerge(cdr w1,w,simplew1)
else if car w1=1 then (car w).sqmerge(cdr w1,cdr w,simplew1) else
append(if simplew1 then list car w1 else car w1,car w).
sqmerge(cdr w1,cdr w,simplew1);
symbolic procedure multup l;
% l is a list of s.f.'s. result is s.f. for product of elements of l;
begin scalar res;
res:=1;
while not null l do <<
res:=multf(res,car l);
l:=cdr l >>;
return res
end;
symbolic procedure diflist(l,cl,x,rl);
% Differentiates l (list of s.f.'s) wrt x to produce the sum of
% terms for the derivative of numr of 1st part of answer. cl is
% coefficient list (s.f.'s) & rl is list of derivatives we have
% dealt with so far. Result is s.q.;
if null l then nil ./ 1
else begin scalar temp;
temp:=!*multf(multup rl,multup cdr l);
temp:=!*multsq(difff(car l,x),!*f2q temp);
temp:=!*multsq(temp,(car cl) ./ 1);
return addsq(temp,diflist(cdr l,cdr cl,x,(car l).rl))
end;
symbolic procedure multsqfree w;
% W is list of sqfree factors. result is product of each LIST IN W
% to give one polynomial for each sqfree power;
if null w then nil
else (multup car w).multsqfree cdr w;
symbolic procedure l2lsf l;
% L is a list of kernels. result is a list of same members as s.f.'s;
if null l then nil
else ((mksp(car l,1) .* 1) .+ nil).l2lsf cdr l;
symbolic procedure dfnumr(x,dl);
% Gives the derivative of the numr of the 1st part of answer.
% dl is list of any exponential or 1+tan**2 that occur in integrand
% denr. these are divided out from result before handing it back.
% result is s.q., ready for printing.
begin scalar temp1,temp2,coeflist,qlist,count;
if not null sqfr then <<
count:=0;
qlist:=cdr sqfr;
coeflist:=nil;
while not null qlist do <<
count:=count+1;
coeflist:=count.coeflist;
qlist:=cdr qlist >>;
coeflist:=reverse coeflist >>;
temp1:=!*multsq(diflist(l2lsf zlist,l2lsf indexlist,x,nil),
!*f2q multup sqfr);
if not null sqfr and not null cdr sqfr then <<
temp2:=!*multsq(diflist(cdr sqfr,coeflist,x,nil),
!*f2q multup l2lsf zlist);
temp2:=!*multsq(temp2,(car sqfr) ./ 1) >>
else temp2:=nil ./ 1;
temp1:=addsq(temp1,negsq temp2);
temp2:=cdr temp1;
temp1:=car temp1;
qlist:=nil;
while not null dl do <<
if not car dl member qlist then qlist:=(car dl).qlist;
dl:=cdr dl >>;
while not null qlist do <<
temp1:=quotf(temp1,car qlist);
qlist:=cdr qlist >>;
return temp1 ./ temp2
end;
symbolic procedure difflogs(ll,denm1,x);
% LL is list of log terms (with coeffts), den is common denominator
% over which they are to be put. Result is s.q. for derivative of all
% these wrt x.
if null ll then nil ./ 1
else begin scalar temp,qu,cvar,logoratan,arg;
logoratan:=caar ll;
cvar:=cadar ll;
arg:=cddar ll;
temp:=!*multsq(cvar ./ 1,diffsq(arg,x));
if logoratan='iden then qu:=1 ./ 1
else if logoratan='log then qu:=arg
else if logoratan='atan then qu:=addsq(1 ./ 1,!*multsq(arg,arg))
else interr "Logoratan=? in difflogs";
%Note call to special division routine;
qu:=fquotf(!*multf(!*multf(denm1,numr temp),
denr qu),numr qu);
%*MUST* GO EXACTLY;
temp:=!*multsq(!*invsq (denr temp ./ 1),qu);
%result of fquotf is a s.q;
return !*addsq(temp,difflogs(cdr ll,denm1,x))
end;
symbolic procedure factorlistlist (w,clogflag);
% W is list of lists of sqfree factors in s.f. result is list of log
% terms required for integral answer. the arguments for each log fn
% are in s.q.;
begin scalar res,x,y;
while not null w do <<
x:=car w;
while not null x do <<
y:=facbypp(car x,varlist);
while not null y do <<
res:=append(int!-fac car y,res);
y:=cdr y >>;
x:=cdr x >>;
w:=cdr w >>;
return res
end;
symbolic procedure facbypp(p,vl);
% Use contents/primitive parts to try to factor p.
if null vl then list p
else begin scalar princilap!-part,co;
co:=contents(p,car vl);
vl:=cdr vl;
if co=1 then return facbypp(p,vl); %this var no help.
princilap!-part:=quotf(p,co); %primitive part.
if princilap!-part=1 then return facbypp(p,vl); %again no help;
return nconc(facbypp(princilap!-part,vl),facbypp(co,vl))
end;
endmodule;
module csolve; % routines to do with the C constants.
% Author: John P. Fitch.
fluid '(ccount cmap cmatrix cval loglist neweqn);
global '(!*trint);
exports backsubst4cs,createcmap,findpivot,printspreadc,printvecsq,
spreadc,subst4eliminateds;
imports nth,interr,!*multf,printsf,printsq,quotf,putv,negf,invsq,
negsq,addsq,multsq,mksp,addf,domainp,pnth;
symbolic procedure findpivot cvec;
% Finds first non-zero element in CVEC and returns its cell number.
% If no such element exists, result is nil.
begin scalar i,x;
i:=1;
x:=getv(cvec,i);
while i<ccount and null x do
<< i:=i+1;
x:=getv(cvec,i) >>;
if null x then return nil;
return i
end;
symbolic procedure subst4eliminatedcs(neweqn,substorder,ceqns);
% Substitutes into NEWEQN for all the C's that have been eliminated so
% far. These are given by CEQNS. SUBSTORDER gives the order of
% substitution as well as the constant multipliers. Result is the
% transformed NEWEQN.
if null substorder then neweqn
else begin scalar nxt,row,cvar,temp;
row:=car ceqns;
nxt:=car substorder;
if null (cvar:=getv(neweqn,nxt)) then
return subst4eliminatedcs(neweqn,cdr substorder,cdr ceqns);
nxt:=getv(row,nxt);
for i:=0 : ccount do
<< temp:=!*multf(nxt,getv(neweqn,i));
temp:=addf(temp,negf !*multf(cvar,getv(row,i)));
putv(neweqn,i,temp) >>;
return subst4eliminatedcs(neweqn,cdr substorder,cdr ceqns)
end;
symbolic procedure backsubst4cs(cs2subst,cs2solve,cmatrix);
% Solves the C-eqns and sets vector CVAL to the C-constant values
% CMATRIX is a list of matrix rows for C-eqns after Gaussian
% elimination has been performed. CS2SOLVE is a list of the remaining
% C's to evaluate and CS2SUBST are the C's we have evaluated already.
if null cmatrix then nil
else begin scalar eqnn,cvar,already,substlist,temp,temp2;
eqnn:=car cmatrix;
cvar:=car cs2solve;
already:=nil ./ 1; % The S.Q. nil.
substlist:=cs2subst;
% Now substitute for previously evaluated c's:
while not null substlist do
<< temp:=car substlist;
if not null getv(eqnn,temp) then
already:=addsq(already,multsq(getv(eqnn,temp) ./ 1,
getv(cval,temp)));
substlist:=cdr substlist >>;
% Now solve for the c given by cvar (any remaining c's assumed zero).
temp:=negsq addsq(getv(eqnn,0) ./ 1,already);
if not null (temp2:=quotf(numr temp,getv(eqnn,cvar))) then
temp:=temp2 ./ denr temp
else temp:=multsq(temp,invsq(getv(eqnn,cvar) ./ 1));
if not null numr temp then putv(cval,cvar,
resimp rootextractsq subs2q temp);
backsubst4cs(reversewoc(cvar . reversewoc cs2subst),
cdr cs2solve,cdr cmatrix)
end;
%**********************************************************************
% Routines to deal with linear equations for the constants C.
%**********************************************************************
symbolic procedure createcmap;
%Sets LOGLIST to list of things of form (LOG C-constant f), where f is
% function linear in one of the z-variables and C-constant is in S.F.
% When creating these C-constant names, the CMAP is also set up and
% returned as the result.
begin scalar i,l,c;
l:=loglist;
i:=1;
while not null l do <<
c:=(int!-gensym1('c) . i) . c;
i:=i+1;
rplacd(car l,((mksp(caar c,1) .* 1) .+ nil) . cdar l);
l:=cdr l >>;
if !*trint then printc ("Constants Map" . c);
return c
end;
symbolic procedure spreadc(eqnn,cvec1,w);
% Sets a vector 'cvec1' to coefficients of c<i> in eqnn.
if domainp eqnn then putv(cvec1,0,addf(getv(cvec1,0),
!*multf(eqnn,w)))
else begin scalar mv,t1,t2;
spreadc(red eqnn,cvec1,w);
mv:=mvar eqnn;
t1:=assoc(mv,cmap); %tests if it is a c var.
if not null t1 then return <<
t1:=cdr t1; %loc in vector for this c.
if not (tdeg lt eqnn=1) then interr "Not linear in c eqn";
t2:=addf(getv(cvec1,t1),!*multf(w,lc eqnn));
putv(cvec1,t1,t2) >>;
t1:=((lpow eqnn) .* 1) .+ nil; %this main var as sf.
spreadc(lc eqnn,cvec1,!*multf(w,t1))
end;
symbolic procedure printspreadc cvec1;
begin
for i:=0 : ccount do <<
prin2 i;
printc ":";
printsf(getv(cvec1,i)) >>;
printc "End of printspreadc output"
end;
%symbolic procedure printvecsq cvec;
%% Print contents of cvec which contains s.q.'s (not s.f.'s).
%% Starts from cell 1 not 0 as above routine (printspreadc).
% begin
% for i:=1 : ccount do <<
% prin2 i;
% printc ":";
% if null getv(cvec,i) then printc "0"
% else printsq(getv(cvec,i)) >>;
% printc "End of printvecsq output"
% end;
endmodule;
module cuberoot; % Cube roots of standard forms.
% Authors: Mary Ann Moore and Arthur C. Norman.
fluid '(cuberootflag);
exports cuberootdf;
imports contentsmv,gcdf,!*multf,nrootn,partialdiff,printdf,quotf,vp2,
mksp,mk!*sq,domainp;
symbolic procedure cuberootsq a;
cuberootf numr a ./ cuberootf denr a;
symbolic procedure cuberootf p;
begin scalar ip,qp;
if null p then return nil;
ip:=cuberootf1 p;
qp:=cdr ip;
ip:=car ip; %respectable and nasty parts of the cuberoot.
if numberp qp and onep qp then return ip; %exact root found.
qp:=list('expt,prepf qp,'(quotient 1 3));
cuberootflag:=t; %symbolic cube-root introduced.
qp:=(mksp(qp,1).* 1) .+ nil;
return !*multf(ip,qp)
end;
symbolic procedure cuberootf1 p;
% Returns a . b with p=a**2*b.
% Does this need power reduction?
if domainp p then nrootn(p,3)
else begin scalar co,ppp,g,pg;
co:=contentsmv(p,mvar p,nil); %contents of p.
ppp:=quotf(p,co); %primitive part.
% now consider ppp=p1*p2**2*p3**3*p4**4*...
co:=cuberootf1(co); %process contents via recursion.
g:=gcdf(ppp,partialdiff(ppp,mvar ppp));
% g=p2*p3**2*p4**3*...
if not domainp g then <<
pg:=quotf(ppp,g);
%pg=p1*p2*p3*p4*...
g:=gcdf(g,partialdiff(g,mvar g));
% g=g3*g4**2*g5**3*...
g:=gcdf(g,pg)>>; %a triple factor of ppp.
if domainp g then pg:=1 . ppp
else <<
pg:=quotf(ppp,!*multf(g,!*multf(g,g))); %what's left.
pg:=cuberootf1 pg; %split that up.
rplaca(pg,!*multf(car pg,g))>>;
%put in the thing found here.
rplaca(pg,!*multf(car pg,car co));
rplacd(pg,!*multf(cdr pg,cdr co));
return pg
end;
endmodule;
module idepend; % Routines for considering dependency among variables.
% Authors: Mary Ann Moore and Arthur C. Norman.
fluid '(taylorvariable);
exports dependspl,dependsp,involvesq,involvsf;
imports taylorp,domainp;
symbolic procedure dependsp(x,v);
if null v then t
else if depends(x,v) then x
else if atom x then if x eq v then x else nil
else if car x = '!*sq then involvesq(cadr x,v)
else if taylorp x
then if v eq taylorvariable then taylorvariable else nil
else begin scalar w;
if x=v then return v;
% Check if a prefix form expression depends on the variable v.
% Note this assumes the form x is in normal prefix notation;
w := x; % preserve the dependency;
x := cdr x; % ready to recursively check arguments;
scan: if null x then return nil; % no dependency found;
if dependsp(car x,v) then return w;
x:=cdr x;
go to scan
end;
symbolic procedure involvesq(sq,term);
involvesf(numr sq,term) or involvesf(denr sq,term);
symbolic procedure involvesf(sf,term);
if domainp sf or null sf then nil
else dependsp(mvar sf,term)
or involvesf(lc sf,term)
or involvesf(red sf,term);
symbolic procedure dependspl(dep!-list,var);
% True if any member of deplist (a list of prefix forms) depends on
% var.
dep!-list
and (dependsp(car dep!-list,var) or dependspl(cdr dep!-list,var));
symbolic procedure taylorp exxpr;
% Sees if a random entity is a taylor expression.
not atom exxpr
and not atom car exxpr
and flagp(taylorfunction exxpr,'taylor);
endmodule;
module df2q; % Conversion from distributive to standard forms.
% Authors: Mary Ann Moore and Arthur C. Norman.
fluid '(indexlist zlist);
exports df2q;
imports addf,gcdf,mksp,!*multf,quotf;
comment We assume that results already have reduced powers, so
that no power substitution is necessary;
symbolic procedure df2q p;
% Converts distributed form P to standard quotient;
begin scalar n,d,gg,w;
if null p then return nil ./ 1;
d:=denr lc p;
w:=red p;
while not null w do <<
gg:=gcdf(d,denr lc w); %get denominator of answer...
d:=!*multf(d,quotf(denr lc w,gg));
%..as lcm of denoms in input
w:=red w >>;
n:=nil; %place to build numerator of answer
while not null p do <<
n:=addf(n,!*multf(xl2f(lpow p,zlist,indexlist),
!*multf(numr lc p,quotf(d,denr lc p))));
p:=red p >>;
return n ./ d
end;
symbolic procedure xl2f(l,z,il);
% L is an exponent list from a D.F., Z is the Z-list,
% IL is the list of indices.
% Value is L converted to standard form. ;
if null z then 1
else if car l=0 then xl2f(cdr l,cdr z,cdr il)
else if not atom car l then
begin scalar temp;
if caar l=0 then temp:= car il
else temp:=list('plus,car il,caar l);
temp:=mksp(list('expt,car z,temp),1);
return !*multf(((temp .* 1) .+ nil),
xl2f(cdr l,cdr z,cdr il))
end
% else if minusp car l then ;
% multsq(invsq (((mksp(car z,-car l) .* 1) .+ nil)), ;
% xl2f(cdr l,cdr z,cdr il)) ;
else !*multf((mksp(car z,car l) .* 1) .+ nil,
xl2f(cdr l,cdr z,cdr il));
endmodule;
module distrib; % Routines for manipulating distributed forms.
% Authors: Mary Ann Moore and Arthur C. Norman.
fluid '(indexlist sqrtlist zlist);
exports dfprintform,multbyarbpowers,negdf,quotdfconst,sub1ind,vp1,
vp2,plusdf,multdf,multdfconst,orddf;
imports interr,addsq,negsq,exptsq,simp,domainp,mk!*sq,addf,
multsq,invsq,minusp,mksp,sub1;
%***********************************************************************
% NOTE: The expressions lt,red,lc,lpow have been used on distributed
% forms as the latter's structure is sufficiently similar to
% s.f.'s. However lc df is a s.q. not a s.f. and lpow df is a
% list of the exponents of the variables. This also makes
% lt df different. Red df is d.f. as expected.
%**********************************************************************;
symbolic procedure plusdf(u,v);
% U and V are D.F.'s. Value is D.F. for 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 numr x then y else (lpow u .* x) .+ y)
(!*addsq(lc u,lc v),plusdf(red u,red v))
else if orddf(lpow u,lpow v) then lt u .+ plusdf(red u,v)
else (lt v) .+ plusdf(u,red v);
symbolic procedure orddf(u,v);
% U and V are the LPOW of a D.F. - i.e. the list of exponents ;
% Value is true if LPOW U '>' LPOW V and false otherwise ;
if null u then if null v then interr "Orddf = case"
else interr "Orddf v longer than u"
else if null v then interr "Orddf u longer than v"
else if exptcompare(car u,car v) then t
else if exptcompare(car v,car u) then nil
else orddf(cdr u,cdr v);
symbolic procedure exptcompare(x,y);
if atom x then if atom y then x>y else nil
else if atom y then t
else car x > car y;
symbolic procedure negdf u;
if null u then nil
else (lpow u .* negsq lc u) .+ negdf red u;
symbolic procedure multdf(u,v);
% U and V are D.F.'s. Value is D.F. for U*V;
% reduces squares of square-roots as it goes;
if null u or null v then nil
else begin scalar y;
%use (a+b)*(c+d) = (a*c) + a*(c+d) + b*(c+d);
y:=multerm(lt u,lt v); %leading terms;
y:=plusdf(y,multdf(red u,v));
y:=plusdf(y,multdf((lt u) .+ nil,red v));
return y
end;
symbolic procedure multerm(u,v);
%multiply two terms to give a D.F.;
begin scalar coef;
coef:=!*multsq(cdr u,cdr v); %coefficient part;
return multdfconst(coef,mulpower(car u,car v))
end;
symbolic procedure mulpower(u,v);
% u and v are exponent lists. multiply corresponding forms;
begin scalar r,s;
r:=addexptsdf(u,v);
if not null sqrtlist then s:=reduceroots(r,zlist);
r:=(r .* (1 ./ 1)) .+ nil;
if not (s=nil) then r:=multdf(r,s);
return r
end;
symbolic procedure reduceroots(r,zl);
begin scalar s;
while not null r do <<
if eqcar(car zl,'sqrt) then
s:=tryreduction(r,car zl,s);
r:=cdr r; zl:=cdr zl >>;
return s
end;
symbolic procedure tryreduction(r,var,s);
begin scalar x;
x:=car r; %current exponent
if not atom x then << r:=x; x:=car r >>; %numeric part
if (x=0) or (x=1) then return s; %no reduction possible
x:=divide(x,2);
rplaca(r,cdr x); %reduce exponent as redorded
x:=car x;
var:=simp cadr var; %sqrt arg as a s q
var:=!*exptsq(var,x);
x:=multdfconst(1 ./ denr var,f2df numr var); %distribute
if s=nil then s:=x
else s:=multdf(s,x);
return s
end;
symbolic procedure addexptsdf(x,y);
% X and Y are LPOW's of D.F. Value is list of sum of exponents;
if null x then if null y then nil else interr "X too long"
else if null y then interr "Y too long"
else exptplus(car x,car y).addexptsdf(cdr x,cdr y);
symbolic procedure exptplus(x,y);
if atom x then if atom y then x+y else list (x+car y)
else if atom y then list (car x +y)
else interr "Bad exponent sum";
symbolic procedure multdfconst(x,u);
% X is S.Q. not involving Z variables of D.F. U. Value is D.F.;
% for X*U;
if (null u) or (null numr x) then nil
else lpow u .* !*multsq(x,lc u) .+ multdfconst(x,red u);
%symbolic procedure quotdfconst(x,u);
% multdfconst(!*invsq x,u);
symbolic procedure f2df p;
% P is standard form. Value is P in D.F.;
if domainp p then dfconst(p ./ 1)
else if mvar p member zlist then
plusdf(multdf(vp2df(mvar p,tdeg lt p,zlist),f2df lc p),
f2df red p)
else plusdf(multdfconst(((lpow p .* 1) .+ nil) ./ 1,f2df lc p),
f2df red p);
symbolic procedure vp1(var,degg,z);
% Takes VAR and finds it in Z (=list), raises it to power DEGG and puts
% the result in exponent list form for use in a distributed form.
if null z then interr "Var not in z-list after all"
else if var=car z then degg.vp2 cdr z
else 0 . vp1(var,degg,cdr z);
symbolic procedure vp2 z;
% Makes exponent list of zeroes.
if null z then nil
else 0 . vp2 cdr z;
symbolic procedure vp2df(var,exprn,z);
% Makes VAR**EXPRN into exponent list and then converts the resulting
% power into a distributed form.
% Special care with square-roots.
if eqcar(var,'sqrt) and (exprn>1) then
mulpower(vp1(var,exprn,z),vp2 z)
else (vp1(var,exprn,z) .* (1 ./ 1)) .+ nil;
symbolic procedure dfconst q;
% Makes a distributed form from standard quotient constant Q;
if numr q=nil then nil
else ((vp2 zlist) .* q) .+ nil;
%df2q moved to a section of its own.
symbolic procedure df2printform p;
%Convert to a standard form good enough for printing.
if null p then nil
else begin
scalar mv,co;
mv:=xl2q(lpow p,zlist,indexlist);
if mv=(1 ./ 1) then <<
co:=lc p;
if denr co=1 then return addf(numr co,
df2printform red p);
co:=mksp(mk!*sq co,1);
return (co .* 1) .+ df2printform red p >>;
co:=lc p;
if not (denr co=1) then mv:=!*multsq(mv,1 ./ denr co);
mv:=mksp(mk!*sq mv,1) .* numr co;
return mv .+ df2printform red p
end;
symbolic procedure xl2q(l,z,il);
% L is an exponent list from a D.F.,Z is the Z-list, IL is the list of
% indices. Value is L converted to standard quotient. ;
if null z then 1 ./ 1
else if car l=0 then xl2q(cdr l,cdr z,cdr il)
else if not atom car l then
begin scalar temp;
if caar l=0 then temp:= car il
else temp:=list('plus,car il,caar l);
temp:=mksp(list('expt,car z,temp),1);
return !*multsq(((temp .* 1) .+ nil) ./ 1,
xl2q(cdr l,cdr z,cdr il))
end
else if minusp car l then
!*multsq(!*invsq(((mksp(car z,-car l) .* 1) .+ nil) ./ 1),
xl2q(cdr l,cdr z,cdr il))
else !*multsq(((mksp(car z,car l) .* 1) .+ nil) ./ 1,
xl2q(cdr l,cdr z,cdr il));
%symbolic procedure sub1ind power;
% if atom power then power-1
% else list sub1 car power;
symbolic procedure multbyarbpowers u;
% Multiplies the ordinary D.F., U, by arbitrary powers
% of the z-variables;
% i-1 j-1 k-1
% i.e. x z z ... so result is D.F. with the exponent list
% 1 2
%appropriately altered to contain list elements instead of numeric ones.
if null u then nil
else ((addarbexptsdf lpow u) .* lc u) .+ multbyarbpowers red u;
symbolic procedure addarbexptsdf x;
% Adds the arbitrary powers to powers in exponent list, X, to produce
% new exponent list. e.g. 3 -> (2) to represent x**3 now becoming :
% 3 i-1 i+2 ;
% x * x = x . ;
if null x then nil
else list exptplus(car x,-1) . addarbexptsdf cdr x;
endmodule;
module divide; % Exact division of standard forms to give a S Q.
% Authors: Mary Ann Moore and Arthur C. Norman.
fluid '(residue sqrtlist zlist);
global '(!*trdiv !*trint);
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 residue;
if (!*trint or !*trdiv) then <<
printc "Dfquotdf called on ";
printdf a; printdf b>>;
a:=dfquotdf1(a,b);
if (!*trint or !*trdiv) then << printc "Quotient given as ";
printdf a >>;
if not null residue then begin
scalar gres,w;
if !*trint or !*trdiv then <<
printc "Residue in dfquotdf =";
printdf residue;
printc "Which should be zero";
w:=residue;
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);
residue:=reversewoc residue;
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;
residue:=xd .+ residue;
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 "Dfquotdf trouble:"; 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;
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>>>>>>;
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 (not zerop car xl)
else involvesp(cdr xl,var,cdr zl);
endmodule;
module driver; % Driving routines for integration program.
% Author: Mary Ann Moore and Arthur C. Norman.
% Modifications by: John P. Fitch.
fluid '(!*backtrace
!*exp
!*gcd
!*keepsqrts
!*mcd
!*nolnr
!*purerisch
!*rationalize
!*sqrt
!*structure
!*uncached
basic!-listofnewsqrts
basic!-listofallsqrts
expression
gaussiani
intvar
listofnewsqrts
listofallsqrts
loglist
sqrt!-intvar
sqrt!-places!-alist
variable
varlist
xlogs
zlist);
global '(!*algint !*failhard !*trint);
exports integratesq,simpint,purge,simpint1;
imports algebraiccase,algfnpl,findzvars,getvariables,interr,printsq,
transcendentalcase,varsinlist,kernp,simpcar,prepsq,mksq,simp,
opmtch,formlnr;
switch algint,nolnr,trint;
% Form is int(expr,var,x1,x2,...);
% meaning is integrate expr wrt var, given that the result may
% contain logs of x1,x2,...
% x1, etc are intended for use when the system has to be helped
% in the case that expr is algebraic.
% Extended arguments x1, x2, etc., are not currently supported.
symbolic procedure simpint u;
% Simplifies an integral. First two components of U are the integrand
% and integration variable respectively. Optional succeeding components
% are log forms for the final integral;
begin scalar ans,expression,variable,loglist,w,
!*purerisch,intvar,listofnewsqrts,listofallsqrts,
sqrtfn,sqrt!-intvar,sqrt!-places!-alist,
basic!-listofallsqrts,basic!-listofnewsqrts;
if atom u or null cdr u then rederr "Not enough arguments for INT";
variable := !*a2k cadr u;
w := cddr u;
if w then rederr "Too many arguments to INT";
listofnewsqrts:= list mvar gaussiani; % Initialize for SIMPSQRT.
listofallsqrts:= list (argof mvar gaussiani . gaussiani);
sqrtfn := get('sqrt,'simpfn);
put('sqrt,'simpfn,'proper!-simpsqrt);
% We need explicit settings of several switches during integral
% evaluation. In addition, the current code cannot handle domains
% like floating point, so we suppress it while the integral is
% calculated. UNCACHED is turned on since integrator does its own
% caching.
begin scalar dmode!*,!*exp,!*gcd,!*keepsqrts,!*mcd,!*sqrt,
!*rationalize,!*structure,!*uncached;
!*keepsqrts := !*sqrt := t;
!*exp := !*gcd := !*mcd := !*structure := !*uncached := t;
dmode!* := nil;
if !*algint
then <<intvar:=variable; % until fix JHD code
% Start a clean slate (in terms of SQRTSAVE) for this integral
sqrt!-intvar:=!*q2f simpsqrti variable;
if (red sqrt!-intvar) or (lc sqrt!-intvar neq 1)
or (ldeg sqrt!-intvar neq 1)
then interr "Sqrt(x) not properly formed"
else sqrt!-intvar:=mvar sqrt!-intvar;
basic!-listofallsqrts:=listofallsqrts;
basic!-listofnewsqrts:=listofnewsqrts;
sqrtsave(basic!-listofallsqrts,basic!-listofnewsqrts,
list(variable . variable))>>;
expression := int!-simp car u;
% loglist := for each x in w collect int!-simp x;
ans := errorset('(integratesq expression variable loglist),
!*backtrace,!*backtrace);
end;
if errorp ans
then return <<put('sqrt,'simpfn,sqrtfn);
if !*failhard then error1();
simpint1(expression . variable . w)>>
else ans := car ans;
expression := sqrtchk numr ans ./ sqrtchk denr ans;
% We now need to check that all simplifications have been done
% but we have to make sure INT is not resimplified.
put('int,'simpfn,'simpiden);
ans := errorset('(resimp expression),t,!*backtrace);
put('int,'simpfn,'simpint);
put('sqrt,'simpfn,sqrtfn);
return if errorp ans then error1() else car ans
end;
symbolic procedure sqrtchk u;
% U is a standard form. Result is another standard form with square
% roots replaced by half powers.
if domainp u then u
else if not eqcar(mvar u,'sqrt)
then addf(multpf(lpow u,sqrtchk lc u),sqrtchk red u)
else addf(multpf(mksp(list('expt,cadr mvar u,'(quotient 1 2)),
ldeg u),
sqrtchk lc u),
sqrtchk red u);
symbolic procedure int!-simp u;
%converts U to canonical form, including the resimplification of
% *sq forms;
subs2 resimp simp!* u;
put('int,'simpfn,'simpint);
symbolic procedure integratesq(integrand,var,xlogs);
begin scalar varlist,zlist;
if !*trint then <<
printc "Integrand is...";
printsq integrand >>;
varlist:=getvariables integrand;
varlist:=varsinlist(xlogs,varlist); %in case more exist in xlogs
zlist:=findzvars(varlist,list var,var,nil); %important kernels
%the next section causes problems with nested exponentials or logs;
begin scalar oldzlist;
while oldzlist neq zlist do <<
oldzlist:=zlist;
foreach zz in oldzlist do
zlist:=findzvars(distexp(pseudodiff(zz,var)),zlist,var,t)>>
end;
if !*trint then <<
printc "with 'new' functions :";
print zlist >>;
if !*purerisch and not allowedfns zlist
then return simpint1 (integrand . var.nil);
% If it is not suitable for Risch;
varlist:=purge(zlist,varlist);
% Now zlist is list of things that depend on x, and varlist is list
% of constant kernels in integrand;
if !*algint and cdr zlist and algfnpl(zlist,var)
then return algebraiccase(integrand,zlist,varlist)
else return transcendentalcase(integrand,var,xlogs,zlist,varlist)
end;
symbolic procedure distexp(l);
if null l then nil
else if atom car l then car l . distexp cdr l
else if (caar l = 'expt) and (cadar l = 'e) then
begin scalar ll;
ll:=caddr car l;
if eqcar(ll,'plus) then <<
ll:=foreach x in cdr ll collect list('expt,'e,x);
return ('times . ll) . distexp cdr l >>
else return car l . distexp cdr l
end
else distexp car l . distexp cdr l;
symbolic procedure pseudodiff(a,var);
if atom a then nil
else if car a memq '(atan equal log plus quotient sqrt times)
then begin scalar aa,bb;
foreach zz in cdr a do <<
bb:=pseudodiff(zz,var);
if aa then aa:=bb . aa else bb >>;
return aa
end
else if car a eq 'expt
then if depends(cadr a,var) then
prepsq simp list('log,cadr a) .
cadr a .
caddr a .
append(pseudodiff(cadr a,var),pseudodiff(caddr a,var))
else caddr a . pseudodiff(caddr a,var)
else list prepsq simpdf(list(a,var));
symbolic procedure simpint1 u;
begin scalar v,!*sqrt;
u := 'int . prepsq car u . cdr u;
if (v := formlnr u) neq u
then if !*nolnr
then <<v:= simp subst('int!*,'int,v);
return remakesf numr v ./ remakesf denr v>>
else <<!*nolnr:= nil . !*nolnr;
u:=errorset(list('simp,mkquote v),
!*backtrace,!*backtrace);
if pairp u then v:=car u;
!*nolnr:= cdr !*nolnr;
return v>>;
return if (v := opmtch u) then simp v
else if !*failhard then rederr "FAILHARD switch set"
else mksq(u,1)
end;
mkop 'int!*;
put('int!*,'simpfn,'simpint!*);
symbolic procedure simpint!* u;
begin scalar x;
return if (x := opmtch('int . u)) then simp x
else simpiden('int!* . u)
end;
symbolic procedure remakesf u;
%remakes standard form U, substituting operator INT for INT!*;
if domainp u then u
else addf(multpf(if eqcar(mvar u,'int!*)
then mksp('int . cdr mvar u,ldeg u)
else lpow u,remakesf lc u),
remakesf red u);
symbolic procedure allowedfns u;
if null u
then t
else if atom car u or
flagp(caar u,'transcendental)
then allowedfns cdr u
else nil;
symbolic procedure purge(a,b);
if null a then b
else if null b then nil
else purge(cdr a,delete(car a,b));
endmodule;
module d3d4; % Splitting of cubics and quartics.
% Authors: Mary Ann Moore and Arthur C. Norman.
fluid '(knowndiscrimsign zlist);
global '(!*trint);
exports cubic,quartic;
imports covecdf,cuberootf,nth,forceazero,makepolydf,multdf,multdfconst,
!*multf,negdf,plusdf,printdf,printsf,quadratic,sqrtf,vp1,vp2,addf,
negf;
symbolic procedure cubic(pol,var,res);
%Split the univariate (wrt z-vars) cubic pol, at least if a
%change of origin puts it in the form (x-a)**3-b=0;
begin scalar a,b,c,d,v,shift,p;
v:=covecdf(pol,var,3);
shift:=forceazero(v,3); %make coeff x**2 vanish.
%also checks univariate.
% if shift='failed then go to prime;
a:=getv(v,3); b:=getv(v,2); %=0, I hope!;
c:=getv(v,1); d:=getv(v,0);
if !*trint then << printc "Cubic has coefficients";
printsf a; printsf b;
printsf c; printsf d >>;
if not null c then <<
if !*trint then printc "Cubic too hard to split";
go to exit >>;
a:=cuberootf(a); %can't ever fail;
d:=cuberootf(d);
if !*trint then << printc "Cube roots of a and d are";
printsf a; printsf d>>;
%now a*(x+shift)+d is a factor of pol;
%create x+shift in p;
p:=(vp2 zlist .* shift) .+ nil;
p:=(vp1(var,1,zlist) .* (1 ./ 1)) .+ p; %(x+shift);
b:=nil;
b:=(vp2 zlist .* (d ./ 1)) .+ b;
b:=plusdf(b,multdfconst(a ./ 1,p));
b:=makepolydf b; %get rid of denominator.
if !*trint then << printc "One factor of the cubic is";
printdf b >>;
res:=('log . b) . res;
%now form the (quadratic) cofactor;
b:=(vp2 zlist .* (!*multf(d,d) ./ 1)) .+ nil;
b:=plusdf(b,multdfconst(negf !*multf(a,d) ./ 1,p));
b:=plusdf(b,multdfconst(!*multf(a,a) ./ 1,
multdf(p,p)));
return quadratic(makepolydf b,var,res); %deal with what is left;
prime:
if !*trint then printc "The following cubic does not split";
exit:
if !*trint then printdf pol;
return ('log . pol) . res
end;
symbolic procedure quartic(pol,var,res);
%Splits univariate (wrt z-vars) quartics that can be written
%in the form (x-a)**4+b*(x-a)**2+c;
begin scalar a,b,c,d,ee,v,shift,p,q,p1,p2,dsc;
v:=covecdf(pol,var,4);
shift:=forceazero(v,4); %make coeff x**3 vanish;
% if shift='failed then go to prime;
a:=getv(v,4); b:=getv(v,3); %=0, I hope.
c:=getv(v,2); d:=getv(v,1);
ee:=getv(v,0);
if !*trint then << printc "Quartic has coefficients";
printsf a; printsf b;
printsf c; printsf d;
printsf ee >>;
if d
then <<if !*trint then printc "Quartic too hard to split";
go to exit >>;
b:=c; c:=ee; %squash up the notation;
if knowndiscrimsign eq 'negative then go to complex;
dsc := addf(!*multf(b,b),multf(-4,!*multf(a,c)));
p2 := minusf c;
if not p2 and minusf dsc then go to complex;
p1 := null b or minusf b;
if not p1 then if p2 then p1 := t else p2 := t;
p1 := if p1 then 'positive else 'negative;
p2 := if p2 then 'negative else 'positive;
a := sqrtf a;
dsc := sqrtf dsc;
if a eq 'failed or dsc eq 'failed then go to prime;
ee := invsq(addf(a,a) ./ 1);
d := multsq(addf(b,negf dsc) ./ 1,ee);
ee := multsq(addf(b,dsc) ./ 1,ee);
if !*trint
then <<printc "Quadratic factors will have coefficients";
printsf a; print 0; printsq d;
printc "or"; printsq ee>>;
p := (vp2 zlist .* shift) .+ nil;
p := (vp1(var,1,zlist) .* (1 ./ 1)) .+ p; %(x+shift);
q := multdf(p,p); %square of same;
q := multdfconst(a ./ 1,q);
p := plusdf(q,(vp2 zlist .* d) .+ nil);
q := plusdf(q,(vp2 zlist .* ee) .+ nil);
if !*trint
then <<printc "Allowing for change of origin:";
printdf p; printdf q>>;
knowndiscrimsign := p1;
res := quadratic(p,var,res);
knowndiscrimsign := p2;
res := quadratic(q,var,res);
go to quarticdone;
complex:
a:=sqrtf(a);
c:=sqrtf(c);
if a eq 'failed or c eq 'failed then go to prime;
b:=addf(!*multf(2,!*multf(a,c)),negf b);
b:=sqrtf b;
if b eq 'failed then go to prime;
%now a*(x+shift)**2 (+/-) b*(x+shift) + c is a factor.
if !*trint
then << printc "Quadratic factors will have coefficients";
printsf a; printsf b; printsf c>>;
p:=(vp2 zlist .* shift) .+ nil;
p:=(vp1(var,1,zlist) .* (1 ./ 1)) .+ p; %(x+shift);
q:=multdf(p,p); %square of same;
p:=multdfconst(b ./ 1,p);
q:=multdfconst(a ./ 1,q);
q:=plusdf(q,(vp2 zlist .* (c ./ 1)) .+ nil);
if !*trint then <<
printc "Allowing for change of origin, p (+/-) q with p,q=";
printdf p; printdf q>>;
%now p+q and p-q are the factors of the quartic;
knowndiscrimsign := 'negative;
res:=quadratic(plusdf(q,p),var,res);
res:=quadratic(plusdf(q,negdf p),var,res);
quarticdone:
knowndiscrimsign := nil;
if !*trint then printc "Quartic done";
return res;
prime:
if !*trint then printc "The following quartic does not split";
exit:
if !*trint then printdf pol;
return ('log . pol) . res
end;
endmodule;
module factr; % Crude factorization routine for integrator.
% Authors: Mary Ann Moore and Arthur C. Norman.
fluid '(zlist);
global '(!*trint);
exports int!-fac,var2df;
imports cubic,df2q,f2df,interr,multdf,printdf,quadratic,quartic,unifac,
uniform,vp1,vp2,sub1;
symbolic procedure int!-fac x;
% Input: primitive, square-free polynomial (s.form).
%output:
% list of 'factors' wrt zlist
% each item in this list is either
% log . sq
% or atan . sq
% and these logs and arctans are all that is needed in the
% integration of 1/(argument).
begin scalar res,pol,dset,var,degree,vars;
pol:=f2df x; %convert to distributed form
dset:=degreeset(pol);
%now extract factors of the form 'x' or 'log(x)' etc;
%these correspond to items in dset with a non-zero cdr.
begin scalar zl,ds;
zl:=zlist; ds:=dset;
while not null ds do <<
if onep cdar ds then <<
res:=('log . var2df(car zl,1,zlist)) . res;
%record in answer.
pol:=multdf(var2df(car zl,-1,zlist),pol);
%divide out.
if !*trint then << printc "Trivial factor found";
printdf cdar res>>;
rplaca(ds,sub1 caar ds . cdar ds) >>
else if null zerop cdar ds then
interr "Repeated trivial factor in arg to factor";
zl:=cdr zl; ds:=cdr ds >>;
end; %single term factors all removed now.
dset:=mapcar(dset,function car); %get lower bounds.
if !*trint
then printc ("Upper bounds of remaining factors are now: " .
dset);
if dset=vp2 zlist then go to finished; %thing left is constant.
begin scalar ds,zl;
var:=car zlist; degree:=car dset;
if not zerop degree then vars:=var . vars;
ds:=cdr dset; zl:=cdr zlist;
while not null ds do <<
if not zerop car ds then <<
vars:=car zl . vars;
if zerop degree or degree>car ds then <<
var:=car zl; degree:=car ds >> >>;
zl:=cdr zl; ds:=cdr ds >>
end;
% Now var is variable that this poly involves to lowest degree.
% Degree is the degree of the poly in same variable.
if !*trint
then printc ("Best var is " . var . "with exponent " .
degree);
if onep degree then <<
res:=('log . pol) . res; %certainly irreducible.
if !*trint
then << printc "The following is certainly irreducible";
printdf pol>>;
go to finished >>;
if degree=2 then <<
if !*trint then << printc "Quadratic";
printdf pol>>;
res:=quadratic(pol,var,res);
go to finished >>;
dset:=uniform(pol,var);
if not (dset='failed) then <<
if !*trint then << printc "Univariate polynomial";
printdf pol >>;
res:=unifac(dset,var,degree,res);
go to finished >>;
if not null cdr vars then go to nasty; %only try univariate now.
if degree=3 then <<
if !*trint then << printc "Cubic";
printdf pol>>;
res:=cubic(pol,var,res);
% if !*overlaymode then excise 'd3d4;
go to finished >>;
if degree=4 then <<
if !*trint then << printc "Quartic";
printdf pol>>;
res:=quartic(pol,var,res);
% if !*overlaymode then excise 'd3d4;
go to finished>>;
%else abandon hope and hand back some rubbish.
nasty:
res:=('log . pol) . res;
printc
"The following polynomial has not been properly factored";
printdf pol;
go to finished;
finished: %res is a list of d.f. s as required
pol:=nil; %convert back to standard forms
while not null res do
begin scalar type,arg;
type:=caar res; arg:=cdar res;
arg:=df2q arg;
if type='log then rplacd(arg,1);
pol:=(type . arg) . pol;
res:=cdr res end;
return pol
end;
symbolic procedure var2df(var,n,zlist);
((vp1(var,n,zlist) .* (1 ./ 1)) .+ nil);
symbolic procedure degreeset pol;
% Finds degree bounds for all vars in distributed form poly.
degreesub(dbl lpow pol,red pol);
symbolic procedure dbl x;
% Converts list of x into list of (x . x).
if null x then nil
else (car x . car x) . dbl cdr x;
symbolic procedure degreesub(cur,pol);
% Update degree bounds 'cur' to include info about pol.
<<
while not null pol do <<
cur:=degreesub1(cur,lpow pol);
pol:=red pol >>;
cur >>;
symbolic procedure degreesub1(cur,nxt);
%Merge information from exponent set next into cur.
if null cur then nil
else degreesub2(car cur,car nxt) . degreesub1(cdr cur,cdr nxt);
symbolic procedure degreesub2(two,one);
max(car two,one) . min(cdr two,one);
endmodule;
module ibasics; % Some basic support routines for integrator.
% Authors: Mary Ann Moore and Arthur C. Norman.
fluid '(!*backtrace !*gcd !*sqfree indexlist sqrtflag sqrtlist
varlist zlist);
global '(!*trint);
exports partialdiff,printdf,rationalintegrate,interr;
imports df2printform,printsf,varsinsf,addsq,multsq,multd,mksp;
symbolic procedure printdf u;
% Print distributed form via cheap conversion to reduce structure.
begin scalar !*gcd;
printsf df2printform u;
end;
symbolic procedure !*n2sq(u1);
if u1=0 then nil . 1 else u1 . 1;
symbolic procedure indx(n);
if n<2 then (list 1) else(n . indx(isub1 n));
symbolic procedure interr mess;
<<(!*trint or !*backtrace)
and <<prin2 "***** INTEGRATION PACKAGE ERROR: "; printc mess>>;
error1()>>;
symbolic procedure rationalintegrate(x,var);
begin scalar n,d;
n:=numr x; d:=denr x;
if not var member varsinsf(d,nil) then
return !*multsq(polynomialintegrate(n,var),1 ./ d);
rederr "Rational integration not coded yet"
end;
symbolic procedure polynomialintegrate(x,v);
% Integrate standard form. result is standard quotient.
if null x then nil ./ 1
else if atom x then ((mksp(v,1) .* 1) .+ nil) ./ 1
else begin scalar r;
r:=polynomialintegrate(red x,v); % deal with reductum
if v=mvar x then begin scalar degree,newlt;
degree:=1+tdeg lt x;
newlt:=((mksp(v,degree) .* lc x) .+ nil) ./ 1; % up exponent
r:=addsq(!*multsq(newlt,1 ./ degree),r)
end
else begin scalar newterm;
newterm:=(((lpow x) .* 1) .+ nil) ./ 1;
newterm:=!*multsq(newterm,polynomialintegrate(lc x,v));
r:=addsq(r,newterm)
end;
return r
end;
symbolic procedure partialdiff(p,v);
% Partial differentiation of p wrt v - p is s.f. as is result.
if domainp p then nil
else
if v=mvar p then
(lambda x; if x=1 then lc p
else ((mksp(v,x-1) .* multd(x,lc p))
.+ partialdiff(red p,v)))
(tdeg lt p)
else
(lambda x; if null x then partialdiff(red p,v)
else ((lpow p .* x) .+ partialdiff(red p,v)))
(partialdiff(lc p,v));
put('pdiff,'simpfn,'simppdiff);
symbolic procedure ncdr(l,n);
% we can use small integer arithmetic here.
if n=0 then l else ncdr(cdr l,isub1 n);
symbolic procedure mkilist(old,term);
if null old then nil
else term.mkilist(cdr old,term);
%symbolic procedure addin(lista,first,listb);
%if null lista
% then nil
% else ((first.car listb).car lista).addin(cdr lista,first,cdr listb);
symbolic procedure removeduplicates(u);
% Purges duplicates from the list passed to it.
if null u then nil
else if (atom u) then u.nil
else if member(car u,cdr u)
then removeduplicates cdr u
else (car u).removeduplicates cdr u;
symbolic procedure jsqfree(sf,var);
begin
varlist:=getvariables(sf ./ 1);
zlist:=findzvars(varlist,list var,var,nil);
sqrtlist:=findsqrts varlist; % before the purge
sqrtflag:=not null sqrtlist;
varlist:=purge(zlist,varlist);
if sf eq !*sqfree
then return list list sf
else return sqfree(sf,zlist)
end;
symbolic procedure jfactor(sf,var);
begin
scalar varlist,zlist,indexlist,sqrtlist,sqrtflag;
scalar prim,answer,u,v; % scalar var2
prim:=jsqfree(sf,var);
indexlist:=createindices zlist;
prim:=factorlistlist (prim,t);
while prim do <<
if caar prim eq 'log then <<
u:=cdar prim;
u:=!*multsq(numr u ./ 1,1 ./ cdr stt(numr u,var));
v:=denr u;
while not atom v do v:=lc v;
if (numberp v) and (0> v)
then u:=(negf numr u) ./ (negf denr u);
answer:=u.answer >>
else if caar prim eq 'atan then <<
% We can ignore this, since we also get LOG (U**2+1) as another term.
>>
else interr "Unexpected term in jfactor";
prim:=cdr prim >>;
return answer
end;
symbolic procedure stt(u,x);
if domainp u
then if u eq nil
then ((-1) . nil)
else (0 . u)
else if mvar u eq x
then ldeg u . lc u
else if ordop(x,mvar u)
then (0 . u)
else begin
scalar ltlc,ltrest;
ltlc:=stt(lc u,x);
ltrest:= stt(red u,x);
if car ltlc = car ltrest then go to merge;
if car ltlc > car ltrest
then return car ltlc .
!*multf(cdr ltlc,(lpow u .* 1) .+ nil)
else return ltrest;
merge:
return car ltlc.addf(cdr ltrest,
!*multf(cdr ltlc,(lpow u .* 1) .+ nil))
end;
symbolic procedure gcdinonevar(u,v,x);
% Gcd of u and v, regarded as polynnmials in x.
if null u
then v
else if null v
then u
else begin
scalar u1,v1,z,w;
u1:=stt(u,x);
v1:=stt(v,x);
loop:
if (car u1) > (car v1)
then go to ok;
z:=u1;u1:=v1;v1:=z;
z:=u;u:=v;v:=z;
ok:
if car v1 iequal 0
then interr "Coprimeness in gcd";
z:=gcdf(cdr u1,cdr v1);
w:=quotf(cdr u1,z);
if (car u1) neq (car v1)
then w:=!*multf(w,!*p2f mksp(x,(car u1)-(car v1)));
u:=addf(!*multf(v,w),
!*multf(u,negf quotf(cdr v1,z)));
if null u
then return v;
u1:=stt(u,x);
go to loop
end;
symbolic procedure mapply(funct,l);
if null l then rederr "Empty list to mapply"
else if null cdr l then car l
else apply(funct,list(car l,mapply(funct,cdr l)));
symbolic procedure !*lcm!*(u,v);
!*multf(u,quotf(v,gcdf(u,v)));
symbolic procedure multsql(u,l);
% Multiplies (!*multsq) each element of l by u.
if null l
then nil
else !*multsq(u,car l).multsql(u,cdr l);
symbolic procedure intersect(x,y);
if null x then nil else if member(car x,y) then
car(x) . intersect(cdr x,y) else
intersect(cdr x,y);
symbolic procedure mapvec(v,f);
begin
scalar n;
n:=upbv v;
for i:=0:n do
apply(f,list getv(v,i))
end;
endmodule;
module jpatches; % Routines for manipulating sf's with power folding.
% Author: James H. Davenport.
fluid '(sqrtflag);
exports !*addsq,!*multsq,!*invsq,!*multf,squashsqrtsq,!*exptsq,!*exptf;
% !*MULTF(A,B) multiplies the polynomials (standard forms) U and V
% in much the same way as MULTF(U,V) would, EXCEPT...
% (1) !*MULTF inhibits the action of OFF EXP and of non-commutative
% multiplications
% (2) Within !*MULTF powers of square roots, and powers of
% exponential kernels are reduced as if substitution rules
% such as FOR ALL X LET SQRT(X)**2=X were being applied;
% Note that !*MULTF comes between MULTF and !*Q2F SUBS2F MULTF in its
% behaviour, and that it is the responsibility of the user to call it
% in sensible places where its services are needed;
% similarly for the other functions defined here;
%symbolic procedure !*addsq(u,v);
%U and V are standard quotients.
% %Value is canonical sum of U and V;
% if null numr u then v
% else if null numr v then u
% else if denr u=1 and denr v=1 then addf(numr u,numr v) ./ 1
% else begin scalar nu,du,nv,dv,x;
% x := gcdf(du:=denr u,dv:=denr v);
% du:=quotf(du,x); dv:=quotf(dv,x);
% nu:=numr u; nv:=numr v;
% u:=addf(!*multf(nu,dv),!*multf(nv,du));
% if u=nil then return nil ./ 1;
% v:=!*multf(du,denr v);
% return !*ff2sq(u,v)
% end;
%symbolic procedure !*multsq(a,b);
% begin
% scalar n,d;
% n:=!*multf(numr a,numr b);
% d:=!*multf(denr a,denr b);
% return !*ff2sq(n,d)
% end;
%symbolic procedure !*ff2sq(a,b);
% begin
% scalar gg;
% if null a then return nil ./ 1;
% gg:=gcdf(a,b);
% if not (gg=1) then <<
% a:=quotf(a,gg);
% b:=quotf(b,gg) >>;
% if minusf b then <<
% a:=negf a;
% b:=negf b >>;
% return a ./ b
% end;
symbolic procedure !*addsq(u,v);
%U and V are standard quotients.
%Value is canonical sum of U and V;
if null numr u then v
else if null numr v then u
else if denr u=1 and denr v=1 then addf(numr u,numr v) ./ 1
else begin scalar du,dv,x,y,z;
x := gcdf(du:=denr u,dv:=denr v);
du:=quotf(du,x); dv:=quotf(dv,x);
y:=addf(!*multf(dv,numr u),!*multf(du,numr v));
if null y then return nil ./ 1;
z:=!*multf(denr u,dv);
if minusf z then <<y := negf y; z := negf z>>;
if x=1 then return y ./ z;
x := gcdf(y,x);
return if x=1 then y ./ z else quotf(y,x) ./ quotf(z,x)
end;
symbolic procedure !*multsq(u,v);
%U and V are standard quotients. Result is the canonical product of
%U and V with surd powers suitably reduced.
if null numr u or null numr v then nil ./ 1
else if denr u=1 and denr v=1 then !*multf(numr u,numr v) ./ 1
else begin scalar w,x,y;
x := gcdf(numr u,denr v);
y := gcdf(numr v,denr u);
w := !*multf(quotf(numr u,x),quotf(numr v,y));
x := !*multf(quotf(denr u,y),quotf(denr v,x));
if minusf x then <<w := negf w; x := negf x>>;
y := gcdf(w,x); % another factor may have been generated.
return if y=1 then w ./ x else quotf(w,y) ./ quotf(x,y)
end;
symbolic procedure !*invsq a;
% Note that several examples (e.g., int(1/(x**8+1),x)) give a more
% compact result when SQRTFLAG is true if SQRT2TOP is not called.
if sqrtflag then sqrt2top invsq a else invsq a;
symbolic procedure !*multf(u,v);
% U and V are standard forms
% Value is SF for U*V;
begin
scalar x,y;
if null u or null v
then return nil
else if u = 1
then return squashsqrt v
else if v = 1
then return squashsqrt u
else if domainp u
then return multd(u,squashsqrt v)
else if domainp v
then return multd(v,squashsqrt u);
x:=mvar u;
y:=mvar v;
if x eq y
then go to c
else if ordop(x,y)
then go to b;
x:=!*multf(u,lc v);
y:=!*multf(u,red v);
return if null x then y
else if not domainp lc v
and mvar u eq mvar lc v
and not atom mvar u
and car mvar u memq '(expt sqrt)
then addf(!*multf(x,!*p2f lpow v),y) % what about noncom?
else makeupsf(lpow v,x,y);
b: x:=!*multf(lc u,v);
y:=!*multf(red u,v);
return if null x then y
else if not domainp lc u
and mvar lc u eq mvar v
and not atom mvar v
and car mvar v memq '(expt sqrt)
then addf(!*multf(!*p2f lpow u,x),y)
else makeupsf(lpow u,x,y);
c: y:=addf(!*multf(list lt u,red v),!*multf(red u,v));
if eqcar(x,'sqrt)
then return addf(squashsqrt y,!*multfsqrt(x,
!*multf(lc u,lc v),ldeg u + ldeg v))
else if eqcar(x,'expt) and prefix!-rational!-numberp caddr x
then return addf(squashsqrt y,!*multfexpt(x,
!*multf(lc u,lc v),ldeg u + ldeg v));
x:=mkspm(x,ldeg u + ldeg v);
return if null x or null (u:=!*multf(lc u,lc v))
then y
else x .* u .+ y
end;
symbolic procedure makeupsf(u,x,y);
% Makes u .* x .+ y except when u is not a valid lpow (because of
% sqrts).
if atom car u or cdr u = 1 then u .* x .+ y
else if caar u eq 'sqrt then addf(!*multfsqrt(car u,x,cdr u),y)
else if <<begin scalar v;
v:=car u;
if car v neq 'expt then return nil;
v:=caddr v;
if atom v then return nil;
return (car v eq 'quotient
and numberp cadr v
and numberp caddr v)
end >>
then addf(!*multfexpt(car u,x,cdr u),y)
else u .* x .+ y;
symbolic procedure !*multfsqrt(x,u,w);
% This code (Due to Norman a& Davenport) squashes SQRT(...)**2.
begin scalar v;
w:=divide(w,2);
v:=!*q2f simp cadr x;
u:=!*multf(u,exptf(v,car w));
if not zerop cdr w then u:=!*multf(u,!*p2f mksp(x,1));
return u
end;
symbolic procedure !*multfexpt(x,u,w);
begin scalar expon,v;
expon:=caddr x;
x:=cadr x;
w:=w * cadr expon;
expon:=caddr expon;
v:=gcdn(w,expon);
w:=w/v;
v:=expon/v;
if not (w > 0) then rederr "Invalid exponent"
else if v = 1
then return !*multf(u,exptf(if numberp x then x
else if atom x then !*k2f x
else !*q2f if car x eq '!*sq
then argof x
else simp x,
w));
expon:=0;
while not (w < v) do <<expon:=expon + 1; w:=w-v>>;
if expon>0 then u:=!*multf(u,exptf(!*q2f simp x,expon));
if w = 0 then return u;
x:=list('expt,x,list('quotient,1,v));
return multf(squashsqrt u,!*p2f mksp(x,w))
end;
symbolic procedure prefix!-rational!-numberp u;
% Tests for m/n in prefix representation.
eqcar(u,'quotient) and numberp cadr u and numberp caddr u;
symbolic procedure squashsqrtsq sq;
!*multsq(squashsqrt numr sq ./ 1,
1 ./ squashsqrt denr sq);
symbolic procedure squashsqrt sf;
if (not sqrtflag) or atom sf or atom mvar sf
then sf
else if car mvar sf eq 'sqrt and ldeg sf > 1
then addf(squashsqrt red sf,!*multfsqrt(mvar sf,lc sf,ldeg sf))
else if car mvar sf eq 'expt
and prefix!-rational!-numberp caddr mvar sf
and ldeg sf >= caddr caddr mvar sf
then addf(squashsqrt red sf,!*multfexpt(mvar sf,lc sf,ldeg sf))
else (lpow sf .* squashsqrt lc sf) .+ squashsqrt red sf;
%remd 'simpx1;
%symbolic procedure simpx1(u,m,n);
% %u,m and n are prefix expressions;
% %value is the standard quotient expression for u**(m/n);
% begin scalar flg,z;
% if null frlis!* or null xn(frlis!*,flatten (m . n))
% then go to a;
% exptp!* := t;
% return !*k2q list('expt,u,if n=1 then m
% else list('quotient,m,n));
% a: if numberp m and fixp m then go to e
% else if atom m then go to b
% else if car m eq 'minus then go to mns
% else if car m eq 'plus then go to pls
% else if car m eq 'times and numberp cadr m and fixp cadr m
% and numberp n
% then go to tms;
% b: z := 1;
% c: if atom u and not numberp u then flag(list u,'used!*);
% u := list('expt,u,if n=1 then m else list('quotient,m,n));
% if not u member exptl!* then exptl!* := u . exptl!*;
% d: return mksq(u,if flg then -z else z); %u is already in lowest
%% %terms;
% e: if numberp n and fixp n then go to int;
% z := m;
% m := 1;
% go to c;
% mns: m := cadr m;
% if !*mcd then return invsq simpx1(u,m,n);
% flg := not flg;
% go to a;
% pls: z := 1 ./ 1;
% pl1: m := cdr m;
% if null m then return z;
% z := multsq(simpexpt list(u,
% list('quotient,if flg then list('minus,car m)
% else car m,n)),
% z);
% go to pl1;
% tms: z := gcdn(n,cadr m);
% n := n/z;
% z := cadr m/z;
% m := retimes cddr m;
% go to c;
% int:z := divide(m,n);
% if cdr z<0 then z:= (car z - 1) . (cdr z+n);
% if 0 = cdr z
% then return simpexpt list(u,car z);
% if n = 2
% then return multsq(simpexpt list(u,car z),
% simpsqrti u);
% return multsq(simpexpt list(u,car z),
% mksq(list('expt,u,list('quotient,1,n)),cdr z))
% end;
symbolic procedure !*exptsq(a,n);
% raise A to the power N using !*MULTSQ;
if n=0 then 1 ./ 1
else if n=1 then a
else if n<0 then !*exptsq(invsq a,-n)
else begin
scalar q,r;
q:=divide(n,2);
r:=cdr q; q:=car q;
q:=!*exptsq(!*multsq(a,a),q);
if r=0 then return q
else return !*multsq(a,q)
end;
symbolic procedure !*exptf(a,n);
% raise A to the power N using !*MULTF;
if n=0 then 1
else if n=1 then a
else begin
scalar q,r;
q:=divide(n,2);
r:=cdr q; q:=car q;
q:=!*exptf(!*multf(a,a),q);
if r=0 then return q
else return !*multf(a,q)
end;
endmodule;
module hacksqrt; % Routines for manipulation of sqrt expressions.
% Author: James H. Davenport.
fluid '(nestedsqrts thisplace);
exports sqrtsintree,sqrtsinsq,sqrtsinsql,sqrtsinsf,sqrtsign;
exports degreenest,sortsqrts;
imports mkvect,interr,getv,dependsp,union;
symbolic procedure sqrtsintree(u,var,slist);
% Adds to slist all the sqrts in the prefix-type tree u.
if atom u
then slist
else if car u eq '!*sq
then union(slist,sqrtsinsq(cadr u,var))
else if car u eq 'sqrt
then if dependsp(argof u,var)
then <<
slist:=sqrtsintree(argof u,var,slist);
% nested square roots
if member(u,slist)
then slist
else u.slist >>
else slist
else sqrtsintree(car u,var,sqrtsintree(cdr u,var,slist));
symbolic procedure sqrtsinsq(u,var);
% Returns list of all sqrts in sq.
sqrtsinsf(denr u,sqrtsinsf(numr u,nil,var),var);
symbolic procedure sqrtsinsql(u,var);
% Returns list of all sqrts in sq list.
if null u
then nil
else sqrtsinsf(denr car u,
sqrtsinsf(numr car u,sqrtsinsql(cdr u,var),var),var);
symbolic procedure sqrtsinsf(u,slist,var);
% Adds to slist all the sqrts in sf.
if domainp u or null u
then slist
else <<
if eqcar(mvar u,'sqrt) and
dependsp(argof mvar u,var) and
not member(mvar u,slist)
then begin
scalar slist2;
slist2:=sqrtsintree(argof mvar u,var,nil);
if slist2
then <<
nestedsqrts:=t;
slist:=union(slist2,slist) >>;
slist:=(mvar u).slist
end;
sqrtsinsf(lc u,sqrtsinsf(red u,slist,var),var) >>;
symbolic procedure easysqrtsign(slist,things);
% This procedure builds a list of all substitutions for all possible
% combinations of square roots in list.
if null slist
then things
else easysqrtsign(cdr slist,
nconc(mapcons(things,(car slist).(car slist)),
mapcons(things,
list(car slist,'minus,car slist))));
symbolic procedure hardsqrtsign(slist,things);
% This procedure fulfils the same role for nested sqrts
% ***assumption: the simpler sqrts come further up the list.
if null slist
then things
else begin
scalar thisplace,answers,pos,neg;
thisplace:=car slist;
answers:=mapcar(things,function (lambda u;sublis(u,thisplace).u));
pos:=mapcar(answers,function (lambda u;(thisplace.car u).(cdr u)));
% pos is sqrt(f) -> sqrt(innersubst f)
neg:=mapcar(answers,
function (lambda u;list(thisplace,'minus,car u).(cdr u)));
% neg is sqrt(f) -> -sqrt(innersubst f)
return hardsqrtsign(cdr slist,nconc(pos,neg))
end;
symbolic procedure degreenest(pf,var);
% Returns the maximum degree of nesting of var
% inside sqrts in the prefix form pf.
if atom pf
then 0
else if car pf eq 'sqrt
then if dependsp(cadr pf,var)
then iadd1 degreenest(cadr pf,var)
else 0
else if car pf eq 'expt
then if dependsp(cadr pf,var)
then if eqcar(caddr pf,'quotient)
then iadd1 degreenest(cadr pf,var)
else degreenest(cadr pf,var)
else 0
else degreenestl(cdr pf,var);
symbolic procedure degreenestl(u,var);
%Returns max degreenest from list of pfs u.
if null u
then 0
else max(degreenest(car u,var),
degreenestl(cdr u,var));
symbolic procedure sortsqrts(u,var);
% Sorts list of sqrts into order required by hardsqrtsign
% (and many other parts of the package).
begin
scalar i,v;
v:=mkvect(10); %should be good enough!
while u do <<
i:=degreenest(car u,var);
if i iequal 0
then interr "Non-dependent sqrt found";
if i > 10
then interr
"Degree of nesting exceeds 10 (recompile with 10 increased)";
putv(v,i,(car u).getv(v,i));
u:=cdr u >>;
u:=getv(v,10);
for i :=9 step -1 until 1 do
u:=nconc(getv(v,i),u);
return u
end;
symbolic procedure sqrtsign(sqrts,x);
if nestedsqrts then hardsqrtsign(sortsqrts(sqrts,x),list nil)
else easysqrtsign(sqrts,list nil);
endmodule;
module kron; % Kronecker factorization of univ. polys over integers.
% Authors: Mary Ann Moore and Arthur C. Norman.
global '(!*trint);
exports linfac,quadfac;
imports zfactor;
% Only linear and quadratic factors are found.
symbolic procedure linfac(w);
trykr(w,'(0 1));
symbolic procedure quadfac(w);
trykr(w,'(-1 0 1));
symbolic procedure trykr(w,points);
%Look for factor of w by evaluation at (points) and use of
% interpolate. Return (fac . cofac) with fac=nil if none
% found and cofac=nil if nothing worthwhile is left.
begin scalar values,attempt;
if null w then return nil . nil;
if (length points > car w) then return w . nil;
%that says if w is already tiny, it is already factored.
values:=mapcar(points,function (lambda x;
evalat(w,x)));
if !*trint then << printc ("At x= " . points);
printc ("p(x)= " . values)>>;
if 0 member values then go to lucky; %(x-1) is a factor!
values:=mapcar(values,function zfactors);
rplacd(values,mapcar(cdr values,function (lambda y;
append(y,mapcar(y,function minus)))));
if !*trint then <<printc "Possible factors go through some of";
print values>>;
attempt:=search4fac(w,values,nil);
if null attempt then attempt:=nil . w;
return attempt;
lucky: %here (x-1) is a factor because p(0) or p(1) or p(-1)
%vanished and cases p(0), p(-1) will have been removed
%elsewhere.
attempt:='(1 1 -1); %the factor
return attempt . testdiv(w,attempt)
end;
symbolic procedure zfactors n;
% Produces a list of all (positive) integer factors of the integer n.
if n=0 then list 0
else if (n:=abs n)=1 then list 1
else combinationtimes zfactor n;
symbolic procedure search4fac(w,values,cv);
% Combinatorial search. cv gets current selected value-set.
% Returns nil if fails, else factor . cofactor.
if null values then tryfactor(w,cv)
else begin scalar ff,q;
ff:=car values; %try all values here
loop: if null ff then return nil; %no factor found
q:=search4fac(w,cdr values,(car ff) . cv);
if null q then << ff:=cdr ff; go to loop>>;
return q
end;
symbolic procedure tryfactor(w,cv);
% Tests if cv represents a factor of w.
begin scalar ff,q;
if null cddr cv then ff:=linethrough(cadr cv,car cv)
else ff:=quadthrough(caddr cv,cadr cv,car cv);
if ff='failed then return nil; %it does not interpolate
q:=testdiv(w,ff);
if q='failed then return nil; %not a factor
return ff . q
end;
symbolic procedure evalat(p,n);
% Evaluate polynomial at integer point n.
begin scalar r;
r:=0;
p:=cdr p;
while not null p do <<
r:=n*r+car p;
p:=cdr p >>;
return r
end;
symbolic procedure testdiv(a,b);
% Quotient a/b or failed.
begin scalar q;
q:=testdiv1(cdr a,car a,cdr b,car b);
if q='failed then return q;
return (car a-car b) . q
end;
symbolic procedure testdiv1(a,da,b,db);
if da<db then begin
check0: if null a then return nil
else if not zerop car a then return 'failed;
a:=cdr a;
go to check0
end
else begin scalar q;
q:=divide(car a,car b);
if zerop cdr q then q:=car q
else return 'failed;
a:=testdiv1(ambq(cdr a,cdr b,q),da-1,b,db);
if a='failed then return a;
return q . a
end;
symbolic procedure ambq(a,b,q);
% A-B*Q with Q an integer.
if null b then a
else ((car a)-(car b)*q) . ambq(cdr a,cdr b,q);
symbolic procedure linethrough(y0,y1);
begin scalar a;
a:=y1-y0;
if zerop a then return 'failed;
if a<0 then <<a:=-a; y0:=-y0 >>;
if onep gcdn(a,y0) then return list(1,a,y0);
return 'failed
end;
symbolic procedure quadthrough(ym1,y0,y1);
begin scalar a,b,c;
a:=divide(ym1+y1,2);
if zerop cdr a then a:=(car a)-y0
else return 'failed;
if zerop a then return 'failed; %linear things already done.
c:=y0;
b:=divide(y1-ym1,2);
if zerop cdr b then b:=car b
else return 'failed;
if not onep gcdn(a,gcd(b,c)) then return 'failed;
if a<0 then <<a:=-a; b:=-b; c:=-c>>;
return list(2,a,b,c)
end;
endmodule;
module lowdeg; % Splitting of low degree polynomials.
% Author: To be determined.
fluid '(clogflag knowndiscrimsign zlist);
global '(!*trint);
exports forceazero,makepolydf,quadratic,covecdf,exponentdf;
imports dfquotdf,gcdf,interr,minusdfp,multdf,multdfconst,!*multf,
negsq,minusp,printsq,multsq,invsq,pnth,nth,mknill,
negdf,plusdf,printdf,printsq,quotf,sqrtdf,var2df,vp2,addsq,sub1;
symbolic procedure covecdf(pol,var,degree);
% Extract coefficients of polynomial wrt var, given a degree-bound
% degree. Result is a lisp vector.
begin scalar v,x,w;
w:=pol;
v:=mkvect(degree);
while not null w do <<
x:=exponentof(var,lpow w,zlist);
if (x<0) or (x>degree) then interr "Bad degree in covecdf";
putv(v,x,lt w . getv(v,x));
w:=red w >>;
for i:=0:degree do putv(v,i,multdf(reversewoc getv(v,i),
var2df(var,-i,zlist)));
return v
end;
symbolic procedure quadratic(pol,var,res);
% Add in to res logs or arctans corresponding to splitting the
% polynomial. Pol given that it is quadratic wrt var.
% Does not assume pol is univariate.
begin scalar a,b,c,w,discrim;
w:=covecdf(pol,var,2);
a:=getv(w,2); b:=getv(w,1); c:=getv(w,0);
% that split the quadratic up to find the coefficients a,b,c.
if !*trint then << printc "a="; printdf a;
printc "b="; printdf b;
printc "c="; printdf c>>;
discrim:=plusdf(multdf(b,b),
multdfconst((-4) . 1,multdf(a,c)));
if !*trint then << printc "Discriminant is";
printdf discrim>>;
if null discrim then interr "Discrim=0 in quadratic";
if knowndiscrimsign
then <<if knowndiscrimsign eq 'negative then go to atancase>>
else if (not clogflag) and (minusdfp discrim)
then go to atancase;
discrim:=sqrtdf(discrim);
if discrim='failed then go to nofactors;
if !*trint then << printc "Square root is";
printdf discrim>>;
w:=var2df(var,1,zlist);
w:=multdf(w,a);
b:=multdfconst(1 ./ 2,b);
discrim:=multdfconst(1 ./ 2,discrim);
w:=plusdf(w,b); %a*x+b/2.
a:=plusdf(w,discrim); b:=plusdf(w,negdf(discrim));
if !*trint then << printc "Factors are";
printdf a; printdf b>>;
return ('log . a) . ('log . b) . res;
atancase:
discrim:=sqrtdf negdf discrim; %sqrt(4*a*c-b**2) this time!
if discrim='failed then go to nofactors; %sqrt did not exist?
res := ('log . pol) . res; %one part of the answer
a:=multdf(a,var2df(var,1,zlist));
a:=plusdf(b,multdfconst(2 ./ 1,a));
a:=dfquotdf(a,discrim); %assumes division is exact
return ('atan . a) . res;
nofactors:
if !*trint
then <<printc
"The following quadratic does not seem to factor";
printdf pol>>;
return ('log . pol) . res
end;
symbolic procedure exponentof(var,l,zl);
if null zl then interr "Var not found in exponentof"
else if var=car zl then car l
else exponentof(var,cdr l,cdr zl);
symbolic procedure df2sf a;
if null a then nil
else if ((null red a) and
(denr lc a = 1) and
(lpow a=vp2 zlist)) then numr lc a
else interr "Nasty cubic or quartic";
symbolic procedure makepolydf p;
% Multiply df by lcm of denominators of all coefficient denominators.
begin scalar h,w;
if null(w:=p) then return nil; %poly is zero already.
h:=denr lc w; %a good start.
w:=red w;
while not null w do <<
h:=quotf(!*multf(h,denr lc w),gcdf(h,denr lc w));
w:=red w >>;
%h is now lcm of denominators.
return multdfconst(h ./ 1,p)
end;
symbolic procedure forceazero(p,n);
%Shift polynomial p so that coeff of x**(n-1) vanishes.
%Return the amount of the shift, update (vector) p.
begin scalar r,i,w;
for i:=0:n do putv(p,i,df2sf getv(p,i)); %convert to polys.
r:=getv(p,n-1);
if null r then return nil ./ 1; %already zero.
r:= !*multsq(r ./ 1,invsq(!*multf(n,getv(p,n)) ./ 1));
% Used to be subs2q multsq
%the shift amount.
%now I have to set p:=subst(x-r,x,p) and then reduce to sf again.
if !*trint then << printc "Shift is by ";
printsq r>>;
w:=mkvect(n); %workspace vector.
for i:=0:n do putv(w,i,nil ./ 1); %zero it.
i:=n;
while not minusp i do <<
mulvecbyxr(w,negsq r,n); %W:=(X-R)*W;
putv(w,0,addsq(getv(w,0),getv(p,i) ./ 1));
i:=i-1 >>;
if !*trint then << printc "SQ shifted poly is";
print w>>;
for i:=0:n do putv(p,i,getv(w,i));
w:=denr getv(p,0);
for i:=1:n do w:=quotf(!*multf(w,denr getv(p,i)),
gcdf(w,denr getv(p,i)));
for i:=0:n do putv(p,i,numr !*multsq(getv(p,i),w ./ 1));
% Used to be subs2q multsq
w:=getv(p,0);
for i:=1:n do w:=gcdf(w,getv(p,i));
if not (w=1) then
for i:=0:n do putv(p,i,quotf(getv(p,i),w));
if !*trint then << printc "Final shifted poly is ";
print p>>;
return r
end;
symbolic procedure mulvecbyxr(w,r,n);
% W is a vector representing a poly of degree n.
% Multiply it by (x+r).
begin scalar i,im1;
i:=n;
im1:=sub1 i;
while not minusp im1 do <<
putv(w,i,!*addsq(getv(w,im1),!*multsq(r,getv(w,i))));
% Used to be subs2q addsq/multsq
i:=im1; im1:=sub1 i >>;
putv(w,0,!*multsq(getv(w,0),r));
% Used to be subs2q multsq
return w
end;
endmodule;
module reform; % Reformulate expressions using C-constant substitution.
% Authors: Mary Ann Moore and Arthur C. Norman.
fluid '(cmap cval loglist ulist);
global '(!*trint);
exports logstosq,substinulist;
imports prepsq,mksp,nth,multsq,addsq,domainp,invsq,plusdf;
symbolic procedure substinulist ulist;
% Substitutes for the C-constants in the values of the U's given in
% ULIST. Result is a D.F.
if null ulist then nil
else begin scalar temp,lcu;
lcu:=lc ulist;
temp:=evaluateuconst numr lcu;
if null numr temp then temp:=nil
else temp:=((lpow ulist) .*
!*multsq(temp,!*invsq(denr lcu ./ 1))) .+ nil;
return plusdf(temp,substinulist red ulist)
end;
symbolic procedure evaluateuconst coefft;
% Substitutes for the C-constants into COEFFT (=S.F.). Result is S.Q.;
if null coefft or domainp coefft then coefft ./ 1
else begin scalar temp;
if null(temp:=assoc(mvar coefft,cmap)) then
temp:=(!*p2f lpow coefft) ./ 1
else temp:=getv(cval,cdr temp);
temp:=!*multsq(temp,evaluateuconst(lc coefft));
% Next line had addsq previously
return !*addsq(temp,evaluateuconst(red coefft))
end;
symbolic procedure logstosq;
% Converts LOGLIST to sum of the log terms as a S.Q.;
begin scalar lglst,logsq,i,temp;
i:=1;
lglst:=loglist;
logsq:=nil ./ 1;
loop: if null lglst then return logsq;
temp:=cddr car lglst;
if !*trint
then <<printc "SF arg for log etc ="; printc temp>>;
if not (caar lglst='iden) then <<
temp:=prepsq temp; %convert to prefix form.
temp:=list(caar lglst,temp); %function name.
temp:=((mksp(temp,1) .* 1) .+ nil) ./ 1 >>;
temp:=!*multsq(temp,getv(cval,i));
% Next line had addsq previously
logsq:=!*addsq(temp,logsq);
lglst:=cdr lglst;
i:=i+1;
go to loop
end;
endmodule;
module simplog; % Simplify logarithms.
% Authors: Mary Ann Moore and Arthur C. Norman.
exports simplog,simplogsq;
imports quotf,prepf,mksp,simp!*,!*multsq,simptimes,addsq,minusf,negf,
addf,comfac,negsq,mk!*sq,carx;
symbolic procedure simplog(exxpr); simplogi carx(exxpr,'simplog);
symbolic procedure simplogi(sq);
begin
if atom sq then go to simplify;
if car sq eq 'times
then return simpplus(for each u in cdr sq
collect mk!*sq simplogi u);
if car sq eq 'quotient
then return addsq(simplogi cadr sq,
negsq simplogi caddr sq);
if car sq eq 'expt
then return simptimes list(caddr sq,
mk!*sq simplogi cadr sq);
if car sq eq 'nthroot
then return !*multsq(1 ./ caddr sq,simplogi cadr sq);
% we had (nthroot of n).
if car sq eq 'sqrt
then return !*multsq(1 ./ 2,simplogi argof sq);
if car sq = '!*sq
then return simplogsq cadr sq;
simplify:
sq:=simp!* sq;
return simplogsq sq
end;
symbolic procedure simplogsq sq;
addsq((simplog2 numr sq),negsq(simplog2 denr sq));
symbolic procedure simplog2(sf);
if atom sf
then if null sf then rederr "Log 0 formed"
else if numberp sf
then if sf iequal 1
then nil ./ 1
else if sf iequal 0 then rederr "Log 0 formed"
else((mksp(list('log,sf),1) .* 1) .+ nil) ./ 1
else formlog(sf)
else begin
scalar form;
form:=comfac sf;
if not null car form
then return addsq(formlog(form .+ nil),
simplog2 quotf(sf,form .+ nil));
% we have killed common powers.
form:=cdr form;
if form neq 1
then return addsq(simplog2 form,
simplog2 quotf(sf,form));
% remove a common factor from the sf.
return (formlog sf)
end;
symbolic procedure formlogterm(sf);
begin
scalar u;
u:=mvar sf;
if not atom u and (car u member '(times sqrt expt nthroot))
then u:=addsq(simplog2 lc sf,
!*multsq(simplogi u,simp!* ldeg sf))
else if (lc sf iequal 1) and (ldeg sf iequal 1)
then u:=((mksp(list('log,u),1) .* 1) .+ nil) ./ 1
else u:=addsq(simptimes list(list('log,u),ldeg sf),
simplog2 lc sf);
return u
end;
symbolic procedure formlog sf;
if null red sf
then formlogterm sf
else if minusf sf
then addf((mksp(list('log,-1),1) .* 1) .+ nil,
formlog2 negf sf) ./ 1
else (formlog2 sf) ./ 1;
symbolic procedure formlog2 sf;
((mksp(list('log,prepf sf),1) .* 1) .+ nil);
endmodule;
module simpsqrt; % Simplify square roots.
% Authors: Mary Ann Moore and Arthur C. Norman.
% Heavily modified J.H. Davenport for algebraic functions.
fluid '(!*backtrace !*conscount !*galois !*pvar basic!-listofallsqrts
gaussiani basic!-listofnewsqrts kord!* knowntobeindep
listofallsqrts listofnewsqrts sqrtflag sqrtlist
sqrt!-places!-alist variable varlist zlist);
global '(!*tra !*trint);
% This module should be rewritten in terms of the REDUCE function
% SIMPSQRT;
% remd 'simpsqrt;
exports proper!-simpsqrt,simpsqrti,simpsqrtsq,simpsqrt2,sqrtsave,
newplace,actualsimpsqrt,formsqrt;
symbolic procedure proper!-simpsqrt(exprn);
simpsqrti carx(exprn,'proper!-simpsqrt);
symbolic procedure simpsqrti sq;
begin
scalar u;
if atom sq
then if numberp sq
then return (simpsqrt2 sq) ./ 1
else if (u:=get(sq,'avalue))
then return simpsqrti cadr u
% BEWARE!!! This is VERY system dependent.
else return simpsqrt2((mksp(sq,1) .* 1) .+ nil) ./ 1;
% If it doesnt have an AVALUE then it is itself;
if car sq eq 'times
then return mapply(function multsq,
mapcar(cdr sq,function simpsqrti));
if car sq eq 'quotient
then return multsq(simpsqrti cadr sq,
invsq simpsqrti caddr sq);
if car sq eq 'expt and numberp caddr sq
then if evenp caddr sq
then return simpexpt list(cadr sq,caddr sq / 2)
else return simpexpt
list(mk!*sq simpsqrti cadr sq,caddr sq);
if car sq = '!*sq
then return simpsqrtsq cadr sq;
return simpsqrtsq tidysqrt simp!* sq
end;
symbolic procedure simpsqrtsq sq;
(simpsqrt2 numr sq) ./ (simpsqrt2 denr sq);
symbolic procedure simpsqrt2 sf;
if minusf sf
then if sf iequal -1
then gaussiani
else begin
scalar u;
u:=negf sf;
if numberp u
then return multf(gaussiani,simpsqrt3 u);
% we cannot negate general expressions for the following reason:
% (%%%thesis remark%%%)
% sqrt(x*x-1) under x->1/x gives sqrt(1-x*x)/x=i*sqrt(x*x-1)/x
% under x->1/x gives x*i*sqrt(-1+1/x*x)=i**2*sqrt(x*x-1)
% hence an abysmal catastrophe;
return simpsqrt3 sf
end
else simpsqrt3 sf;
symbolic procedure simpsqrt3 sf;
begin
scalar u;
u:=assoc(sf,listofallsqrts);
if u
then return cdr u;
% now see if 'knowntobeindep'can help.
u:=atsoc(listofnewsqrts,knowntobeindep);
if null u
then go to no;
u:=assoc(sf,cdr u);
if u
then <<
listofallsqrts:=u.listofallsqrts;
return cdr u >>;
no:
u:=actualsimpsqrt sf;
listofallsqrts:=(sf.u).listofallsqrts;
return u
end;
symbolic procedure sqrtsave(u,v,place);
begin
scalar a;
%u is new value of listofallsqrts, v of new.
a:=assoc(place,sqrt!-places!-alist);
if null a
then sqrt!-places!-alist:=(place.(listofnewsqrts.listofallsqrts))
.sqrt!-places!-alist
else rplacd(a,listofnewsqrts.listofallsqrts);
listofnewsqrts:=v;
% throw away things we are not going to need in future.
if not !*galois
then listofallsqrts:=u;
% we cannot guarantee the validity of our calculations.
if listofallsqrts eq u
then return nil;
v:=listofallsqrts;
while not (cdr v eq u) do
v:=cdr v;
rplacd(v,nil);
% listofallsqrts is all those added since routine was entered.
v:=atsoc(listofnewsqrts,knowntobeindep);
if null v
then knowntobeindep:=(listofnewsqrts.listofallsqrts)
. knowntobeindep
else rplacd(v,union(cdr v,listofallsqrts));
listofallsqrts:=u;
return nil
end;
symbolic procedure newplace(u);
% Says to restart algebraic bases at a new place u.
begin
scalar v;
v:=assoc(u,sqrt!-places!-alist);
if null v
then <<
listofallsqrts:=basic!-listofallsqrts;
listofnewsqrts:=basic!-listofnewsqrts >>
else <<
v:=cdr v;
listofnewsqrts:=car v;
listofallsqrts:=cdr v >>;
return if v then v
else listofnewsqrts.listofallsqrts
end;
symbolic procedure mknewsqrt u;
% U is prefix form.
begin
scalar v,w;
if not !*galois then go to new;
% no checking required.
v:=addf(!*p2f mksp(!*pvar,2),negf !*q2f tidysqrt simp u);
% count !*conscount;
w:=errorset(list('afactor,mkquote v,mkquote !*pvar),t,!*backtrace);
% if !*tra then <<
% princ "*** That took ";
% princ (!*conscount - count nil);
% printc " conses" >>;
% count 0;
if atom w
then go to new
else w:=car w; %the actual result of afactor.
if cdr w
then go to notnew;
new:
w:=sqrtify u;
listofnewsqrts:=w . listofnewsqrts;
return !*kk2f w;
notnew:
w:=car w;
v:=stt(w,!*pvar);
if car v neq 1
then rederr "HELP ...";
w:=addf(w,multf(cdr v,(mksp(!*pvar,car v) .* -1) .+nil));
w:=sqrt2top(w ./ cdr v);
w:=quotf(numr w,denr w);
if null w
then rederr "Division failure";
return w
end;
symbolic procedure actualsimpsqrt(sf);
if sf iequal -1
then gaussiani
else actualsqrtinner(sf,listofnewsqrts);
symbolic procedure actualsqrtinner(sf,l);
if null l
then actualsimpsqrt2 sf
else begin
scalar z;
% z:=simp argof car l;
% if denr z neq 1 or (z := numr z) iequal -1
z:=!*q2f simp argof car l;
if z iequal -1
then return actualsqrtinner(sf,cdr l);
z:=quotf(sf,z);
if null z
then return actualsqrtinner(sf,cdr l);
return !*multf(!*kk2f car l,actualsimpsqrt z)
end;
symbolic procedure actualsimpsqrt2(sf);
if atom sf
then if null sf
then nil
else if numberp sf
then if sf < 0
then multf(gaussiani,actualsimpsqrt2(- sf))
%Above 2 lines inserted JHD 4 Sept 80;
% test case: SQRT(B*X**2-C)/SQRT(X);
else begin
scalar n;
n:=int!-sqrt sf;
% Changed for conformity with DEC20 LISP JHD July 1982;
if not fixp n
then return mknewsqrt sf
else return n
end
else mknewsqrt(sf)
else begin
scalar form;
form:=comfac sf;
if car form
then return multf((if null cdr sf and (car sf = form)
then formsqrt(form .+ nil)
else simpsqrt2(form .+ nil)),
%The above 2 lines changed by JHD;
%(following suggestions of Morrison);
%to conform to Standard LISP 4 Sept 80;
simpsqrt2 quotf(sf,form .+ nil));
% we have killed common powers.
form:=cdr form;
if form neq 1
then return multf(simpsqrt2 form,
simpsqrt2 quotf(sf,form));
% remove a common factor from the sf.
return formsqrt sf
end;
symbolic procedure int!-sqrt n;
% Return sqrt of n if same is exact, or something non-numeric
% otherwise.
if not numberp n then 'nonnumeric
else if n<0 then 'negative
else if floatp n then sqrt!-float n
else if n<2 then n
else int!-nr(n,(n+1)/2);
symbolic procedure int!-nr(n,root);
% root is an overestimate here. nr moves downwards to root;
begin
scalar w;
w:=root*root;
if n=w then return root;
w:=(root+n/root)/2;
if w>=root then return !*q2f simpsqrt list n;
return int!-nr(n,w)
end;
symbolic procedure formsqrt(sf);
if (null red sf)
then if (lc sf iequal 1) and (ldeg sf iequal 1)
then mknewsqrt mvar sf
else multf(if evenp ldeg sf
then !*p2f mksp(mvar sf,ldeg sf / 2)
else exptf(mknewsqrt mvar sf,ldeg sf),simpsqrt2 lc sf)
else begin
scalar varlist,zlist,sqrtlist,sqrtflag;
scalar v,l,n,w;
% This returns a list, the i-th member of which is
% a list of the factors of multiplicity i (as s.f's);
v:=jsqfree(sf,if variable and involvesf(sf,variable)
then variable
else findatom mvar sf);
% VARIABLE is the best thing to do square-free
% decompositions with respect to, but anything
% else will do if VARIABLE is not set;
if null cdr v and null cdar v
then return mknewsqrt prepf sf;
% The JSQFREE did nothing;
l:=nil;
n:=0;
while v do <<
n:=n+1;
w:=car v;
while w do <<
l:=list('expt,mk!*sq !*f2q car w,n) . l;
w:=cdr w >>;
v:=cdr v >>;
if null cdr l
then l:=car l
else l:='times.l;
% makes L into a valid prefix form;
return !*q2f simpsqrti l
end;
symbolic procedure findatom pf;
if atom pf
then pf
else findatom argof pf;
symbolic procedure sqrtify u;
% Actually creates the sqrt.
begin
scalar v,n,prinlist;
n:=degreenest(u,nil);
% v:=list('sqrt,u);
v := mksqrt u; % To ensure sqrt**2 rule set.
prinlist:=member(v,kord!*);
if prinlist then return car prinlist;
% This might improve sharing.
% anyway, we mustn't re-add this object to KORD!*;
while kord!* and
eqcar(car kord!*,'sqrt) and
(degreenest(argof car kord!*,nil) > n) do <<
prinlist:=(car kord!*) . prinlist;
kord!*:=cdr kord!* >>;
kord!*:=v . kord!*;
while prinlist do <<
kord!*:=(car prinlist) . kord!*;
prinlist:=cdr prinlist >>;
return v
end;
% deflist ('((sqrt (((x) quotient (sqrt x) (times 2 x))))),'dfn);
% put('sqrt,'simpfn,'proper!-simpsqrt); % Now in driver.
endmodule;
module isolve; % Routines for solving the final reduction equation.
% Author: Mary Ann Moore and Arthur C. Norman.
% Modifications by: John P. Fitch.
fluid '(badpart
ccount
cmap
cmatrix
cval
indexlist
lhs!*
lorder
nnn
orderofelim
pt
rhs!*
sillieslist
tanlist
ulist
zlist);
global '(!*number!* !*statistics !*trint);
exports solve!-for!-u;
imports nth,findpivot,gcdf,int!-gensym1,mkvect,interr,multdfconst,
!*multf!*,negdf,orddf,plusdf,printdf,printsf,printspreadc,printsq,
quotf,putv,spreadc,subst4eliminatedcs,mknill,pnth,domainp,addf,
invsq,multsq;
symbolic procedure uterm(powu,rhs!*);
% Finds the contribution from RHS!* of reduction equation, of the;
% U-coefficient given by POWU. Result is in D.F.;
if null rhs!* then nil
else begin scalar coef,power;
power:=addinds(powu,lpow rhs!*);
coef:=evaluatecoeffts(numr lc rhs!*,powu);
if null coef then return uterm(powu,red rhs!*);
coef:=coef ./ denr lc rhs!*;
return plusdf((power .* coef) .+ nil,uterm(powu,red rhs!*))
end;
symbolic procedure solve!-for!-u(rhs!*,lhs!*,ulist);
% Solves the reduction eqn LHS!*=RHS!*. Returns list of U-coefficients
% and their values (ULIST are those we have so far), and a list of
% C-equations to be solved (CLIST are the eqns we have so far).
if null lhs!* then ulist
else begin scalar u,lpowlhs;
lpowlhs:=lpow lhs!*;
begin scalar ll,mm,chge; ll:=maxorder(rhs!*,zlist,0);
mm:=lorder;
while mm do << if car ll < car mm then
<< chge:=t; rplaca(mm,car ll) >>;
ll:=cdr ll; mm:=cdr mm >>;
if !*trint and chge then << print ("Maxorder now ".lorder) >>
end;
u:=pickupu(rhs!*,lpow lhs!*,t);
if null u then
<< if !*trint then << printc "***** C-equation to solve:";
printsf numr lc lhs!*;
printc " = 0";
printc " ">>;
% Remove a zero constant from the lhs, rather than use
% Gauss Elim;
if gausselimn(numr lc lhs!*,lt lhs!*) then
lhs!*:=squashconstants(red lhs!*)
else lhs!*:=red lhs!* >>
else
<< ulist:=(car u .
!*multsq(coefdf(lhs!*,lpowlhs),invsq cdr u)).ulist;
% used to be subs2q multsq
if !*statistics then !*number!*:=!*number!*+1;
if !*trint then <<prin2 "***** U"; prin2 car u; prin2t " =";
printsq multsq(coefdf(lhs!*,lpowlhs),invsq cdr u);
printc " ">>;
lhs!*:=plusdf(lhs!*,
negdf multdfconst(cdar ulist,uterm(car u,rhs!*))) >>;
if !*trint then << printc ".... LHS is now:";
printdf lhs!*;
printc " ">>;
return solve!-for!-u(rhs!*,lhs!*,ulist)
end;
symbolic procedure squashconstants(express);
begin scalar constlst,ii,xp,cl,subby,cmt,xx;
constlst:=reverse cmap;
cmt:=cmatrix;
xxx: xx:=car cmt; % Look at next row of Cmatrix;
cl:=constlst; % and list of the names;
ii:=1; % will become index of removed constant;
while not getv(xx,ii) do
<< ii:=ii+1; cl:=cdr cl >>;
subby:=caar cl; %II is now index, and SUBBY the name;
if member(subby,sillieslist) then
<<cmt:=cdr cmt; go to xxx>>; %This loop must terminate;
% This is because at least one constant remains;
xp:=prepsq !*f2q getv(xx,0); % start to build up the answer;
cl:=cdr cl;
if not (ccount=ii) then for jj:=ii+1:ccount do <<
if getv(xx,jj) then
xp:=list('plus,xp,
list('times,caar cl,
prepsq !*f2q getv(xx,jj)));
cl:=cdr cl >>;
xp:=list('quotient,list('minus,xp),
prepsq !*f2q getv(xx,ii));
if !*trint then << prin2 "Replace "; prin2 subby;
prin2 " by "; printsq simp xp >>;
sillieslist:=subby . sillieslist;
return subdf(express,xp,subby)
end;
symbolic procedure checku(ulist,u);
% Checks that U is not already in ULIST - ie. that this u-coefficient;
% has not already been given a value;
if null ulist then nil
else if (car u) = caar ulist then t
else checku(cdr ulist,u);
symbolic procedure checku1(powu,rhs!*);
%Checks that use of a particular U-term will not cause trouble;
%by introducing negative exponents into lhs when it is used;
begin
top:
if null rhs!* then return nil;
if negind(powu,lpow rhs!*) then
if not null evaluatecoeffts(numr lc rhs!*,powu) then return t;
rhs!*:=red rhs!*;
go to top
end;
symbolic procedure negind(pu,pr);
%check if substituting index values in power gives rise to -ve
% exponents;
if null pu then nil
else if (car pu+caar pr)<0 then t
else negind(cdr pu,cdr pr);
symbolic procedure evaluatecoeffts(coefft,indlist);
% Substitutes the values of the i,j,k,...'s that appear in the S.F. ;
% COEFFT (=coefficient of r.h.s. of reduction equation). Result is S.F.;
if null coefft or domainp coefft then
if coefft=0 then nil else coefft
else begin scalar temp;
if mvar coefft member indexlist then
temp:=valuecoefft(mvar coefft,indlist,indexlist)
else temp:=!*p2f lpow coefft;
temp:=!*multf(temp,evaluatecoeffts(lc coefft,indlist));
return addf(temp,evaluatecoeffts(red coefft,indlist))
end;
symbolic procedure valuecoefft(var,indvalues,indlist);
% Finds the value of VAR, which should be in INDLIST, given INDVALUES;
% - the corresponding values of INDLIST variables;
if null indlist then interr "Valuecoefft - no value"
else if var eq car indlist then
if car indvalues=0 then nil
else car indvalues
else valuecoefft(var,cdr indvalues,cdr indlist);
symbolic procedure addinds(powu,powrhs);
% Adds indices in POWU to those in POWRHS. Result is LPOW of D.F.;
if null powu then if null powrhs then nil
else interr "Powrhs too long"
else if null powrhs then interr "Powu too long"
else (car powu + caar powrhs).addinds(cdr powu,cdr powrhs);
symbolic procedure pickupu(rhs!*,powlhs,flg);
% Picks up the 'lowest' U coefficient from RHS!* if it exists and
% returns it in the form of LT of D.F..
% Returns NIL if no legal term in RHS!* can be found.
% POWLHS is the power we want to match (LPOW of D.F).
% and COEFFU is the list of previous coefficients that must be zero;
begin scalar coeffu,u;
pt:=rhs!*;
top:
if null pt then return nil; %no term found - failed;
u:=nextu(lt pt,powlhs); %check this term...;
if null u then go to notthisone;
if not testord(car u,lorder) then go to neverthisone;
if not checkcoeffts(coeffu,car u) then go to notthisone;
%that inhibited clobbering things already passed over;
if checku(ulist,u) then go to notthisone;
%that avoided redefining a u value;
if checku1(car u,rhs!*) then go to neverthisone;
%avoid introduction of negative exponents;
if flg then
u:=patchuptan(list u,powlhs,red pt,rhs!*);
return u;
neverthisone:
coeffu:=(lc pt) . coeffu;
notthisone:
pt:=red pt;
go to top
end;
symbolic procedure patchuptan(u,powlhs,rpt,rhs!*);
begin
scalar uu,cc,dd,tanlist,redu,redu1;
pt:=rpt;
while pt do <<
if (uu:=pickupu(pt,powlhs,nil))
and testord(car uu,lorder) then <<
% Nasty found, patch it up;
cc:=(int!-gensym1('!C).caar u).cc;
% CC is an alist of constants;
if !*trint then <<prin2 "***** U";
prin2 caar u;
prin2t " =";
print caar cc >>;
redu:=plusdf(redu,
multdfconst(!*k2q caar cc,uterm(caar u,rhs!*)));
u:=uu.u
>>;
if pt then pt:=red pt >>;
redu1:=redu;
while redu1 do begin scalar xx; xx:=car redu1;
if !*trint then << prin2 "Introduced residue "; print xx >>;
if (not testord(car xx,lorder)) then <<
if !*trint then <<
printsq cdr xx; printc " = 0" >>;
if dd:=killsingles(cadr xx,cc) then <<
redu:=subdf(redu,0,car dd);
redu1:=subdf(redu1,0,car dd);
ulist:=((cdr dd).(nil ./ 1)).ulist;
u:=rmve(u,cdr dd);
cc:=purgeconst(cc,dd) >>
else redu1:=cdr redu1 >>
else redu1:=cdr redu1 end;
for each xx in redu do <<
if (not testord(car xx,lorder)) then <<
while cc do <<
addctomap(caar cc);
ulist:=((cdar cc).(!*k2q caar cc))
. ulist;
if !*statistics
then !*number!*:=!*number!*+1;
cc:=cdr cc >>;
gausselimn(numr lc redu,lt redu)>> >>;
if redu then << while cc do << addctomap(caar cc);
ulist:=((cdar cc).(!*k2q caar cc)).ulist;
if !*statistics then !*number!*:=!*number!*+1;
cc:=cdr cc >>;
lhs!*:=plusdf(lhs!*,negdf redu) >>;
return car u
end;
symbolic procedure killsingles(xx,cc);
if atom xx then nil
else if not (cdr xx eq nil) then nil
else begin scalar dd;
dd:=assoc(caaar xx,cc);
if dd then return dd;
return killsingles(cdar xx,cc)
end;
symbolic procedure rmve(l,x);
if caar l=x then cdr l else cons(car l,rmve(cdr l,x));
symbolic procedure subdf(a,b,c);
% SUBSTITUTE B FOR C INTO THE DF A;
% Used to get rid of silly constants introduced;
if a=nil then nil else
begin scalar x;
x:=subs2q subf(numr lc a,list (c . b)) ;
if x=(nil . 1) then return subdf(red a,b,c)
else return plusdf(
list ((lpow a).((car x).!*multf(cdr x,denr lc a))),
subdf(red a,b,c))
end;
symbolic procedure testord(a,b);
% Test order of two DF's in recursive fashion;
if null a then t
else if car a leq car b then testord(cdr a,cdr b)
else nil;
symbolic procedure tanfrom(rhs!*,z,nnn);
% We notice that in all bad cases we have (j-num)tan**j...;
% Extract the num;
begin scalar n,zz,r,rr;
r:=rhs!*;
n:=0; zz:=zlist;
while car zz neq z do << n:=n+1; zz:=cdr zz >>;
a: if null r then go to b;
rr:=caar r; % The list of powers;
for i:=1:n do rr:=cdr rr;
if fixp caar rr then if caar rr>0 then <<
rr:=numr cdar r;
if null red rr then rr:=nil ./ 1
else if fixp (rr:=quotf(red rr,lc rr))
then rr:=-rr else rr:=0>>;
if atom rr then go to b;
r:=cdr r;
go to a;
b: if null r then return maxfrom(lhs!*,nnn)+1;
return max(rr,maxfrom(lhs!*,nnn)+1)
end;
symbolic procedure coefdf(y,u);
if y=nil then nil
else if lpow y=u then lc y
else coefdf(red y,u);
symbolic procedure purgeconst(a,b);
% Remove a const from and expression. May be the same as DELETE?;
if null a then nil
else if car a=b then purgeconst(cdr a,b)
else cons(car a,purgeconst(cdr a,b));
symbolic procedure maxorder(rhs!*,z,n);
% Find a limit on the order of terms, theis is ad hoc;
if null z then nil
else if eqcar(car z,'sqrt) then
cons(1,maxorder(rhs!*,cdr z,n+1))
else if (atom car z) or (caar z neq 'tan) then
cons(maxfrom(lhs!*,n)+1,maxorder(rhs!*,cdr z,n+1))
else cons(tanfrom(rhs!*,car z,n),maxorder(rhs!*,cdr z,n+1));
symbolic procedure maxfrom(l,n);
% Largest order in the nth varable;
if null l then 0
else max(nth(caar l,n+1),maxfrom(cdr l,n));
symbolic procedure copy u;
if atom u then u
else cons(copy car u,copy cdr u);
symbolic procedure addctomap cc;
begin
scalar ncval;
ccount:=ccount+1;
ncval:=mkvect(ccount);
for i:=0:(ccount-1) do putv(ncval,i,getv(cval,i));
putv(ncval,ccount,nil ./ 1);
cval:=ncval;
cmap:=(cc . ccount).cmap;
if !*trint then << prin2 "Constant map changed to "; print cmap >>;
cmatrix:=mapcar(cmatrix,function addtovector);
end;
symbolic procedure addtovector v;
begin scalar vv;
vv:=mkvect(ccount);
for i:=0:(ccount-1) do putv(vv,i,getv(v,i));
putv(vv,ccount,nil);
return vv
end;
symbolic procedure checkcoeffts(cl,indv);
% checks to see that the coefficients in CL (coefficient list - S.Q.s);
% are zero when the i,j,k,... are given values in INDV (LPOW of;
% D.F.). if so the result is true else NIL=false;
if null cl then t
else begin scalar res;
res:=evaluatecoeffts(numr car cl,indv);
if not(null res or res=0) then return nil
else return checkcoeffts(cdr cl,indv)
end;
symbolic procedure nextu(ltrhs,powlhs);
% picks out the appropriate U coefficients for term: LTRHS to match the
% powers of the z-variables given in POWLHS (= exponent list of D.F.).
% return this coefficient in form LT of D.F. If U coefficient does
% not exist then result is NIL. If it is multiplied by a zero then
% result is NIL;
if null ltrhs then nil
else begin scalar indlist,ucoefft;
indlist:=subtractinds(powlhs,car ltrhs,nil);
if null indlist then return nil;
ucoefft:=evaluatecoeffts(numr cdr ltrhs,indlist);
if null ucoefft or ucoefft=0 then return nil;
return indlist .* (ucoefft ./ denr cdr ltrhs)
end;
symbolic procedure subtractinds(powlhs,l,sofar);
% subtract the indices in list L from those in POWLHS to find;
% appropriate values for i,j,k,... when equating coefficients of terms;
% on lhs of reduction eqn. SOFAR is the resulting value list we;
% have constructed so far. if any i,j,k,... value is -ve then result;
% is NIL;
if null l then reversewoc sofar
else if ((car powlhs)-(caar l))<0 then nil
else subtractinds(cdr powlhs,cdr l,
((car powlhs)-(caar l)) . sofar);
symbolic procedure gausselimn(equation,tokill);
% Performs Gaussian elimination on the matrix for the c-equations;
% as each c-equation is found. EQUATION is the next one to deal with;
begin scalar newrow,pivot;
if zerop ccount then go to noway; %failure
newrow:=mkvect(ccount);
spreadc(equation,newrow,1);
subst4eliminatedcs(newrow,reverse orderofelim,reverse cmatrix);
pivot:=findpivot newrow;
if null pivot then go to nopivotfound;
orderofelim:=pivot . orderofelim;
newrow:=makeprim newrow; %remove hcf from new equation
cmatrix:=newrow . cmatrix;
% if !*trint then printspreadc newrow;
return t;
nopivotfound:
if null getv(newrow,0) then <<
if !*trint then printc "Already included";
return nil>>; %equation was 0=0
noway:
badpart:=tokill . badpart; %non-integrable term.
if !*trint then printc "Inconsistent";
return nil
end;
symbolic procedure makeprim row;
begin scalar g;
g:=getv(row,0);
for i:=1:ccount do g:=gcdf(g,getv(row,i));
if g neq 1 then
for i:=0:ccount do putv(row,i,quotf(getv(row,i),g));
for i := 0:ccount do
<<g := getv(row,i);
if g and not domainp g
then putv(row,i,numr resimp((rootextractf g) ./ 1))>>;
return row
end;
endmodule;
module sqrtf; % Square root of standard forms.
% Authors: Mary Ann Moore and Arthur C. Norman.
fluid '(!*noextend zlist);
exports minusdfp,sqrtdf,nrootn,domainp,minusf;
imports contentsmv,gcdf,interr,!*multf,partialdiff,printdf,quotf,
simpsqrt2,vp2;
symbolic procedure minusdfp a;
% Test sign of leading coedd of d.f.
if null a then interr "Minusdfp 0 illegal"
else minusf numr lc a;
symbolic procedure sqrtdf l;
% Takes square root of distributive form. "Failed" usually means
% that the square root is not among already existing objects.
if null l then nil
else begin scalar c;
if lpow l=vp2 zlist then go to ok;
c:=sqrtsq df2q l;
if numr c eq 'failed
then return 'failed;
if denr c eq 'failed
then return 'failed;
return for each u in f2df numr c
collect (car u).!*multsq(cdr u,1 ./ denr c);
ok: c:=sqrtsq lc l;
if numr c eq 'failed or
denr c eq 'failed
then return 'failed
else return (lpow l .* c) .+nil
end;
symbolic procedure sqrtsq a;
sqrtf numr a ./ sqrtf denr a;
symbolic procedure sqrtf p;
begin scalar ip,qp;
if null p then return nil;
ip:=sqrtf1 p;
qp:=cdr ip;
ip:=car ip; %respectable and nasty parts of the sqrt.
if qp=1 then return ip; %exact root found.
if !*noextend then return 'failed;
% We cannot add new square roots in this case, since it is
% then impossible to determine if one square root depends
% on another if new ones are being added all the time.
if zlistp qp then return 'failed;
% Liouville's theorem tells you that you never need to add
% new algebraics depending on the variable of integration.
qp:=simpsqrt2 qp;
return !*multf(ip,qp)
end;
symbolic procedure zlistp qp;
if atom qp then member(qp,zlist)
else or(member(mvar qp,zlist),zlistp lc qp,zlistp red qp);
symbolic procedure sqrtf1 p;
% Returns a . b with p=a**2*b.
if domainp p
then if fixp p then nrootn(p,2)
else !*q2f simpsqrt list prepf p . 1
else begin scalar co,pp,g,pg;
co:=contentsmv(p,mvar p,nil); %contents of p.
pp:=quotf(p,co); %primitive part.
co:=sqrtf1(co); %process contents via recursion.
g:=gcdf(pp,partialdiff(pp,mvar pp));
pg:=quotf(pp,g);
g:=gcdf(g,pg); %a repeated factor of pp.
if g=1 then pg:=1 . pp
else <<
pg:= quotf(pp,!*multf(g,g)); %what is still left.
pg:=sqrtf1(pg); %split that up.
rplaca(pg,!*multf(car pg,g))>>;
%put in the thing found here.
rplaca(pg, !*multf(car pg,car co));
rplacd(pg, !*multf(cdr pg,cdr co));
return pg
end;
% NROOTN removed as in REDUCE base.
endmodule;
module tdiff; % Differentiation routines for integrator.
% Authors: Mary Ann Moore and Arthur C. Norman.
exports !-!-simpdf;
imports simpcar,kernp,diffsq,prepsq,msgpri;
flag('(!-!-simpdf),'lose);
symbolic procedure !-!-simpdf u;
% U is a list of forms, the first an expression and the remainder
% kernels and numbers.
% Value is derivative of first form wrt rest of list.
begin scalar v,x,y,tt;
tt := time(); %start the clock;
v := cdr u;
u := simpcar u;
a: if null v or null numr u then go to exit;
x := if null y or y=0 then simpcar v else y;
if null kernp x then go to e;
x := caaaar x;
v := cdr v;
if null v then go to c;
y := simpcar v;
if null numr y then go to d
else if not denr y=1 or not numberp numr y then go to c;
y := car y;
v := cdr v;
b: if y=0 then go to a;
u := diffsq(u,x);
y := y-1;
go to b;
c: u := diffsq(u,x);
go to a;
d: y := nil;
v := cdr v;
go to a;
exit:
print list('time,time()-tt);
return u;
e: msgpri("Differentiation wrt",prepsq x,"not allowed",nil,t)
end;
put('tdf,'simpfn,'!-!-simpdf);
endmodule;
module tidysqrt; % General tidying up of square roots.
% Authors: Mary Ann Moore and Arthur C. Norman.
% Modifications by J.H. Davenport.
exports sqrt2top,tidysqrt;
%symbolic procedure tidysqrtdf a;
% if null a then nil
% else begin scalar tt,r;
% tt:=tidysqrt lc a;
% r:=tidysqrtdf red a;
% if null numr tt then return r;
% return ((lpow a) .* tt) .+ r
% end;
%
symbolic procedure tidysqrt q;
begin scalar nn,dd;
nn:=tidysqrtf numr q;
if null nn then return nil ./ 1; %answer is zero.
dd:=tidysqrtf denr q;
return multsq(nn,invsq dd)
end;
symbolic procedure tidysqrtf p;
%Input - standard form.
%Output - standard quotient.
%Simplifies sqrt(a)**n with n>1.
if domainp p then p ./ 1
else begin scalar v,w;
v:=lpow p;
if car v='i then v:=mksp('(sqrt -1),cdr v); %I->sqrt(-1);
if eqcar(car v,'sqrt) and not onep cdr v
then begin scalar x;
%here we have a reduction to apply.
x:=divide(cdr v,2); %halve exponent.
w:=exptsq(simp cadar v,car x); %rational part of answer.
if not zerop cdr x then w:=multsq(w,
((mksp(car v,1) .* 1) .+ nil) ./ 1);
%the next line allows for the horrors of nested sqrts.
w:=tidysqrt w
end
else w:=((v .* 1) .+ nil) ./ 1;
v:=multsq(w,tidysqrtf lc p);
return addsq(v,tidysqrtf red p)
end;
symbolic procedure multoutdenr q;
% Move sqrts in a sq to the numerator.
begin scalar n,d,root,conj;
n:=numr q;
d:=denr q;
while (root:=findsquareroot d) do <<
conj:=conjugatewrt(d,root);
n:=!*multf(n,conj);
d:=!*multf(d,conj) >>;
while (root:=findnthroot d) do <<
conj:=conjugateexpt(d,root,kord!*);
n:=!*multf(n,conj);
d:=!*multf(d,conj) >>;
return (n . d);
end;
symbolic procedure conjugateexpt(d,root,kord!*);
begin scalar ord,ans,repl,xi;
ord:=caddr caddr root; % the denominator of the exponent;
ans:=1;
kord!*:= (xi:=gensym()) . kord!*;
% XI is an ORD'th root of unity;
for i:=1:ord-1 do <<
ans:=!*multf(ans,numr subf(d,
list(root . list('times,root,list('explt,xi,i)))));
while (mvar ans eq xi) and ldeg ans > ord do
ans:=addf(red ans,(xi) to (ldeg ans - ord) .* lc ans .+ nil);
if (mvar ans eq xi) and ldeg ans = ord then
ans:=addf(red ans,lc ans) >>;
if (mvar ans eq xi) and ldeg ans = ord-1 then <<
repl:=-1;
for i:=1:ord-2 do
repl:=(xi) to i .* -1 .+ repl;
ans:=addf(red ans,!*multf(lc ans,repl)) >>;
if not domainp ans and mvar ans eq xi
then interr "Conjugation failure";
return ans;
end;
symbolic procedure sqrt2top q;
begin
scalar n,d;
n:=multoutdenr q;
d:=denr n;
n:=numr n;
if d eq denr q
then return q;%no change.
if d iequal 1
then return (n ./ 1);
q:=gcdcoeffsofsqrts n;
if q iequal 1
then if minusf d
then return (negf n ./ negf d)
else return (n ./ d);
q:=gcdf(q,d);
n:=quotf(n,q);
d:=quotf(d,q);
if minusf d
then return (negf n ./ negf d)
else return (n ./ d)
end;
%symbolic procedure denrsqrt2top q;
%begin
% scalar n,d;
% n:=multoutdenr q;
% d:=denr n;
% n:=numr n;
% if d eq denr q
% then return d; % no changes;
% if d iequal 1
% then return 1;
% q:=gcdcoeffsofsqrts n;
% if q iequal 1
% then return d;
% q:=gcdf(q,d);
% if q iequal 1
% then return d
% else return quotf(d,q)
% end;
symbolic procedure findsquareroot p;
% Locate a sqrt symbol in poly p.
if domainp p then nil
else begin scalar w;
w:=mvar p; %check main var first.
if atom w
then return nil; %we have passed all sqrts.
if eqcar(w,'sqrt) then return w;
w:=findsquareroot lc p;
if null w then w:=findsquareroot red p;
return w
end;
symbolic procedure findnthroot p;
nil; % Until corrected.
symbolic procedure x!-findnthroot p;
% Locate an n-th root symbol in poly p.
if domainp p then nil
else begin scalar w;
w:=mvar p; %check main var first.
if atom w
then return nil; %we have passed all sqrts.
if eqcar(w,'expt) and eqcar(caddr w,'quotient) then return w;
w:=findnthroot lc p;
if null w then w:=findnthroot red p;
return w
end;
symbolic procedure conjugatewrt(p,var);
% Var -> -var in form p.
if domainp p then p
else if mvar p=var then begin
scalar x,c,r;
x:=tdeg lt p; %degree
c:=lc p; %coefficient
r:=red p; %reductum
x:=remainder(x,2); %now just 0 or 1.
if x=1 then c:=negf c; %-coefficient.
return (lpow p .* c) .+ conjugatewrt(r,var) end
else if ordop(var,mvar p) then p
else (lpow p .* conjugatewrt(lc p,var)) .+
conjugatewrt(red p,var);
symbolic procedure gcdcoeffsofsqrts u;
if atom u
then if numberp u and minusp u
then -u
else u
else if eqcar(mvar u,'sqrt)
then begin
scalar v;
v:=gcdcoeffsofsqrts lc u;
if v iequal 1
then return v
else return gcdf(v,gcdcoeffsofsqrts red u)
end
else begin
scalar root;
root:=findsquareroot u;
if null root
then return u;
u:=makemainvar(u,root);
root:=gcdcoeffsofsqrts lc u;
if root iequal 1
then return 1
else return gcdf(root,gcdcoeffsofsqrts red u)
end;
endmodule;
module trcase; % Driving routine for integration of transcendental fns.
% Authors: Mary Ann Moore and Arthur C. Norman.
% Modifications by: John P. Fitch.
fluid '(!*backtrace
!*nowarnings
!*purerisch
!*reverse
badpart
ccount
cmap
cmatrix
content
cuberootflag
cval
denbad
denominator
indexlist
lhs!*
loglist
lorder
orderofelim
rhs!*
sillieslist
sqfr
sqrtflag
sqrtlist
tanlist
svar
varlist
xlogs
zlist);
% !*reverse: flag to re-order zlist.
% !*nowarnings: flag to lose messages.
global '(!*failhard
!*number!*
!*ratintspecial
!*seplogs
!*spsize!*
!*statistics
!*trint
gensymcount);
switch failhard;
exports transcendentalcase;
imports backsubst4cs,countz,createcmap,createindices,df2q,dfnumr,
difflogs,fsdf,factorlistlist,findsqrts,findtrialdivs,gcdf,mkvect,
interr,logstosq,mergin,multbyarbpowers,!*multf,multsqfree,
printdf,printsq,quotf,rationalintegrate,putv,
simpint1,solve!-for!-u,sqfree,sqmerge,sqrt2top,substinulist,trialdiv,
mergein,negsq,addsq,f2df,mknill,pnth,invsq,multsq,domainp,mk!*sq,
mksp,prettyprint,prepsq;
% Note that SEPLOGS keeps logarithmic part of result together as a
% kernel form, but this can lead to quite messy results.
symbolic
procedure transcendentalcase(integrand,svar,xlogs,zlist,varlist);
begin scalar divlist,jhd!-content,content,prim,sqfr,dfu,indexlist,
% JHD!-CONTENT is local, while CONTENT is free (set in SQFREE).
sillieslist,originalorder,wrongway,
sqrtlist,tanlist,loglist,dflogs,eprim,dfun,unintegrand,
sqrtflag,badpart,rhs!*,lhs!*,gcdq,cmap,cval,orderofelim,cmatrix;
scalar cuberootflag,ccount,denominator,result,denbad;
gensymcount:=0;
integrand:=sqrt2top integrand; % Move the sqrts to the numerator.
if !*trint then << printc "Extension variables z<i> are";
print zlist>>;
if !*ratintspecial and null cdr zlist then
return rationalintegrate(integrand,svar);
% *** now unnormalize integrand, maybe ***.
begin scalar w,gg;
gg:=1;
foreach z in zlist do <<
w:=subs2 diffsq(simp z,svar);
gg:=!*multf(gg,quotf(denr w,gcdf(denr w,gg))) >>;
gg:=quotf(gg,gcdf(gg,denr integrand));
unintegrand:=(!*multf(gg,numr integrand)
./ !*multf(gg,denr integrand));
if !*trint then <<
printc "Unnormalized integrand =";
printsq unintegrand >> end;
divlist:=findtrialdivs zlist;
%also puts some things on loglist sometimes.
% if !*trint
% then << printc "Exponentials and tans to try dividing:";
% print divlist>>;
sqrtlist:=findsqrts zlist;
% if !*trint then << printc "Square-root z-variables";
% print sqrtlist >>;
divlist:=trialdiv(denr unintegrand,divlist);
% if !*trint then << printc "Divisors:";
% print car divlist;
% print cdr divlist>>;
%n.b. the next line also sets 'content' as a free variable.
% Since SQFREE may be used later, we copy it into JHD!-CONTENT.
prim:=sqfree(cdr divlist,zlist);
jhd!-content:=content;
printfactors(prim,nil);
eprim:=sqmerge(countz car divlist,prim,nil);
printfactors(eprim,t);
% if !*trint then << terpri();
% printsf denominator;
% terpri();
% printc "...content is:";
% printsf jhd!-content>>;
sqfr:=for each u in eprim collect multup u;
% sqfr:=multsqfree eprim;
% if !*trint then << printc "...sqfr is:";
% superprint sqfr>>;
if !*reverse then zlist:=reverse zlist; %ALTER ORDER FUNCTION;
indexlist:=createindices zlist;
% if !*trint then << printc "...indices are:";
% superprint indexlist>>;
dfu:=dfnumr(svar,car divlist);
% if !*trint then << terpri();
% printc "************ Derivative of u is:";
% printsq dfu>>;
loglist:=append(loglist,factorlistlist (prim,nil));
loglist:=mergein(xlogs,loglist);
loglist:=mergein(tanlist,loglist);
cmap:=createcmap();
ccount:=length cmap;
if !*trint then << printc "Loglist ";
print loglist >>;
dflogs:=difflogs(loglist,denr unintegrand,svar);
if !*trint then << printc "************ 'Derivative' of logs is:";
printsq dflogs>>;
dflogs:=addsq((numr unintegrand) ./ 1,negsq dflogs);
% Put everything in reduction eqn over common denominator.
gcdq:=gcdf(denr dflogs,denr dfu);
dfun:= !*multf(numr dfu,
denbad:=quotf(denr dflogs,gcdq));
denbad:=!*multf(denr dfu,denbad);
denbad:= !*multf(denr unintegrand,denbad);
dflogs:= !*multf(numr dflogs,quotf(denr dfu,gcdq));
dfu:=dfun;
% Now DFU and DFLOGS are S.F.s.
rhs!*:=multbyarbpowers f2df dfu;
if checkdffail(rhs!*,svar) then interr "Simplification failure";
if !*trint then << printc "Distributed Form of U is:";
printdf rhs!*>>;
lhs!*:=f2df dflogs;
if checkdffail(lhs!*,svar) then interr "Simplification failure";
if !*trint then << printc "Distributed Form of l.h.s. is:";
printdf lhs!*;
terpri()>>;
cval:=mkvect(ccount);
for i:=0 : ccount do putv(cval,i,nil ./ 1);
lorder:=maxorder(rhs!*,zlist,0);
originalorder:=lorder;
if !*trint then << printc "Maximum order determined as ";
print lorder >>;
if !*statistics then << !*number!*:=0;
!*spsize!*:=1;
foreach xx in lorder do
!*spsize!*:=!*spsize!* * (xx+1) >>;
% That calculates the largest U that can appear.
dfun:=solve!-for!-u(rhs!*,lhs!*,nil);
backsubst4cs(nil,orderofelim,cmatrix);
% if !*trint then if not (ccount=0) then printvecsq cval;
if !*statistics then << prin2 !*number!*; prin2 " used out of ";
printc !*spsize!* >>;
badpart:=substinulist badpart;
%substitute for c<i> still in badpart.
dfun:=df2q substinulist dfun;
% if !*trint then superprint dfun;
result:= !*multsq(dfun,!*invsq(denominator ./ 1));
result:= !*multsq(result,!*invsq(jhd!-content ./ 1));
% if !*trint then superprint result;
dflogs:=logstosq();
if not null numr dflogs then <<
if !*seplogs and (not domainp numr result) then <<
result:=mk!*sq result;
result:=(mksp(result,1) .* 1) .+ nil;
result:=result ./ 1 >>;
result:=addsq(result,dflogs)>>;
if !*trint then << superprint result;
terpri();
printc
"*****************************************************";
printc
"************ THE INTEGRAL IS : **********************";
printc
"*****************************************************";
terpri();
printsq result;
terpri()>>;
if not null badpart then <<
if !*trint then printc "plus a bad part";
lhs!*:=badpart;
lorder:=maxorder(rhs!*,zlist,0);
while lorder do <<
if car lorder > car originalorder then
wrongway:=t;
lorder:=cdr lorder;
originalorder:=cdr originalorder >>;
dfun:=df2q badpart;
if !*trint
then <<printsq dfun; printc "Denbad = "; printsf denbad>>;
dfun:= !*multsq(dfun,invsq(denbad ./ 1));
if wrongway then << result:= nil ./ 1; dfun:=integrand >>;
if rootcheckp(unintegrand,svar) then
return simpint1(integrand . svar.nil)
else if !*purerisch or allowedfns zlist then
dfun:=simpint1 (dfun . svar.nil)
else << !*purerisch:=t;
if !*trint
then <<printc " [Transforming ..."; printsq dfun>>;
denbad:=transform(dfun,svar);
if denbad=dfun
then dfun:=simpint1(dfun . svar.nil)
else <<denbad:=errorset('(integratesq denbad svar xlogs),
!*backtrace,!*backtrace);
if not atom denbad then dfun:=untan car denbad
else dfun:=simpint1(dfun . svar.nil) >> >>;
if !*trint then printsq dfun;
if !*failhard then rederr "FAILHARD switch set";
if !*seplogs and not domainp result then <<
result:=mk!*sq result;
if not numberp result
then result:=(mksp(result,1) .* 1) .+ nil;
result:=result ./ 1>>;
result:=addsq(result,dfun) >>;
% if !*overlaymode then excise transcode;
return sqrt2top result
end;
symbolic procedure checkdffail(u,v);
u and (depends(lc u,v) or checkdffail(red u,v));
symbolic procedure printfactors(w,prdenom);
% W is a list of factors to each power. If PRDENOM is true
% this prints denominator of answer, else prints square-free
% decomposition.
begin scalar i,wx;
i:=1;
if prdenom then <<
denominator:=1;
if !*trint
then printc "Denominator of 1st part of answer is:";
if not null w then w:=cdr w >>;
loopx: if w=nil then return;
if !*trint then <<prin2 "Factors of multiplicity "; print i>>;
wx:=car w;
while not null wx do <<
if !*trint then printsf car wx;
for j:=1 : i do
denominator:= !*multf(car wx,denominator);
wx:=cdr wx >>;
i:=i+1;
w:=cdr w;
go to loopx
end;
% unfluid '(dfun svar xlogs);
endmodule;
module halfangle; % Routines for conversion to half angle tangents.
% Author: Steve Harrington.
% Modifications by: John P. Fitch.
fluid '(!*gcd);
exports halfangle,untan;
symbolic procedure transform(u,x);
% Transform the SQ U to remove the 'bad' functions sin, cos, cot etc
% in favor of half angles;
halfangle(u,x);
symbolic procedure quotqq(u1,v1);
multsq(u1, invsq(v1));
symbolic procedure !*subtrq(u1,v1);
addsq(u1, negsq(v1));
symbolic procedure !*int2qm(u1);
if u1=0 then nil . 1 else u1 . 1;
symbolic procedure halfangle(r,x);
% Top level procedure for converting;
% R is a rational expression to be converted,
% X the integration variable.
% A rational expression is returned.
quotqq(hfaglf(numr(r),x), hfaglf(denr(r),x));
symbolic procedure hfaglf(p,x);
% Converting polynomials, a rational expression is returned.
if domainp(p) then !*f2q(p)
else subs2q addsq(multsq(exptsq(hfaglk(mvar(p),x), ldeg(p)),
hfaglf(lc(p),x)),
hfaglf(red(p),x));
symbolic procedure hfaglk(k,x);
% Converting kernels, a rational expression is returned.
begin
scalar kt;
if atom k or not member(x,flatten(cdr(k))) then return !*k2q k;
k := car(k) . hfaglargs(cdr(k), x);
kt := simp list('tan, list('quotient, cadr(k), 2));
return if car(k) = 'sin
then quotqq(multsq(!*int2qm(2),kt), addsq(!*int2qm(1),
exptsq(kt,2)))
else if car(k) = 'cos
then quotqq(!*subtrq(!*int2qm(1),exptsq(kt,2)),addsq(!*int2qm(1),
exptsq(kt,2)))
else if car(k) = 'tan
then quotqq(multsq(!*int2qm(2),kt), !*subtrq(!*int2qm(1),
exptsq(kt,2)))
else if car(k) = 'sinh then
quotqq(!*subtrq(exptsq(!*k2q('expt.('e. cdr k)),2),
!*int2qm(1)), multsq(!*int2qm(2), !*k2q('expt . ('e . cdr(k)))))
else if car(k) = 'cosh then
quotqq(addsq(exptsq(!*k2q('expt.('e. cdr k)),2),
!*int2qm(1)), multsq(!*int2qm(2), !*k2q('expt . ('e . cdr(k)))))
else if car(k) = 'tanh then
quotqq(!*subtrq(exptsq(!*k2q('expt.('e. cdr k)),2),
!*int2qm(1)), addsq(exptsq(!*k2q ('expt.('e.cdr(k))),2),
!*int2qm(1)))
else !*k2q(k); % additional transformation might be added here.
end;
symbolic procedure hfaglargs(l,x);
% Conversion of argument list.
if null l then nil
else prepsq(hfaglk(car(l),x)) . hfaglargs(cdr(l), x);
symbolic procedure untanf x;
% This should be done by a table.
begin scalar y,z,w;
if domainp x then return x . 1;
y := mvar x;
if eqcar(y,'int) then error1(); % assume all is hopeless.
z := ldeg x;
w := 1 . 1;
y :=
if atom y then !*k2q y
else if car y eq 'tan
then if evenp z
then <<z := z/2;
simp list('quotient,
list('plus,
list('minus,
list('cos,
'times
. (2 . cdr y))),
1),list('plus,
list('cos,
'times
. (2 . cdr y)),
1))>>
else if z=1
then simp list('quotient,
list('plus,
list('minus,
list('cos,
'times . (2 . cdr y))),
1),list('sin,
'times . (2 . cdr y)))
else <<z := (z - 1)/2;
w :=
simp list('quotient,
list('plus,
list('minus,
list('cos,
'times
. (2 . cdr y))),
1),list('sin,
'times
. (2 . cdr y)));
simp list('quotient,
list('plus,
list('minus,
list('cos,
'times
. (2 . cdr y))),
1),list('plus,
list('cos,
'times
. (2 . cdr y)),
1))>>
else simp y;
return addsq(multsq(multsq(exptsq(y,z),untanf lc x),w),
untanf red x)
end;
symbolic procedure untanlist(y);
if null y then nil
else (prepsq (untan(simp car y)) . untanlist(cdr y));
symbolic procedure untan(x);
% Expects x to be canonical quotient.
begin scalar y;
y:=cossqchk sinsqrdchk multsq(untanf(numr x),
invsq untanf(denr x));
return if length flatten y>length flatten x then x else y
end;
symbolic procedure sinsqrdchk(x);
multsq(sinsqchkf(numr x), invsq sinsqchkf(denr x));
symbolic procedure sinsqchkf(x);
begin
scalar y,z,w;
if domainp x then return x . 1;
y := mvar x;
z := ldeg x;
w := 1 . 1;
y := if eqcar(y,'sin) then if evenp z
then <<z := quotient(z,2);
simp list('plus,1,list('minus,
list('expt,('cos . cdr(y)),2)))>>
else if z = 1 then !*k2q y
else << z := quotient(difference(z,1),2); w := !*k2q y;
simp list('plus,1,list('minus,
list('expt,('cos . cdr(y)),2)))>>
else !*k2q y;
return addsq(multsq(multsq(exptsq(y,z),sinsqchkf(lc x)),w),
sinsqchkf(red x));
end;
symbolic procedure cossqchkf(x);
begin
scalar y,z,w,x1,x2;
if domainp x then return x . 1;
y := mvar x;
z := ldeg x;
w := 1 . 1;
x1 := cossqchkf(lc x);
x2 := cossqchkf(red x);
x := addsq(multsq(!*p2q lpow x,x1),x2);
y := if eqcar(y,'cos) then if evenp z
then <<z := quotient(z,2);
simp list('plus,1,list('minus,
list('expt,('sin . cdr(y)),2)))>>
else if z = 1 then !*k2q y
else << z := quotient(difference(z,1),2); w := !*k2q y;
simp list('plus,1,list('minus,
list('expt,('sin . cdr(y)),2)))>>
else !*k2q y;
y := addsq(multsq(multsq(exptsq(y,z),w),x1),x2);
return if length(y) > length(x) then x else y;
end;
symbolic procedure cossqchk(x);
begin scalar !*gcd;
!*gcd := t;
return multsq(cossqchkf(numr x), invsq cossqchkf(denr x))
end;
symbolic procedure lrootchk(l,x);
% Checks each member of list l for a root.
if null l then nil else krootchk(car l, x) or lrootchk(cdr l, x);
symbolic procedure krootchk(f,x);
% Checks a kernel to see if it is a root.
if atom f then nil
else if car(f) = 'sqrt and member(x, flatten cdr f) then t
else if car(f) = 'expt
and not atom caddr(f)
and caaddr(f) = 'quotient
and member(x, flatten cadr f) then t
else lrootchk(cdr f, x);
symbolic procedure rootchk1p(f,x);
% Checks polynomial for a root.
if domainp f then nil
else krootchk(mvar f,x) or rootchk1p(lc f,x) or rootchk1p(red f,x);
symbolic procedure rootcheckp(f,x);
% Checks rational (standard quotient) for a root.
rootchk1p(numr f,x) or rootchk1p(denr f,x);
endmodule;
module trialdiv; % Trial division routines.
% Authors: Mary Ann Moore and Arthur C. Norman.
fluid '(denominator loglist tanlist);
global '(!*trint);
exports countz,findsqrts,findtrialdivs,trialdiv,simp,mksp;
imports !*multf,printsf,quotf;
symbolic procedure countz dl;
% DL is a list of S.F.s;
begin scalar s,n,rl;
loop2: if null dl then return arrangelistz rl;
n:=1;
loop1: n:=n+1;
s:=car dl;
dl:=cdr dl;
if not null dl and (s eq car dl) then
go to loop1
else rl:=(s.n).rl;
go to loop2
end;
symbolic procedure arrangelistz d;
begin scalar n,s,rl,r;
n:=1;
if null d then return rl;
loopd: if (cdar d)=n then s:=(caar d).s
else r:=(car d).r;
d:=cdr d;
if not null d then go to loopd;
d:=r;
rl:=s.rl;
s:=nil;
r:=nil;
n:=n+1;
if not null d then go to loopd;
return reversewoc rl
end;
symbolic procedure findtrialdivs zl;
%zl is list of kernels found in integrand. result is a list
%giving things to be treated specially in the integration
%viz: exps and tans.
%Result is list of form ((a . b) ...)
% with a a kernel and car a=expt or tan
% and b a standard form for either expt or (1+tan**2).
begin scalar dlists1,args1;
while not null zl do <<
if exportan car zl then <<
if caar zl='tan
then << args1:=(mksp(car zl,2) .* 1) .+ 1;
tanlist:=(args1 ./ 1) . tanlist>>
else args1:=!*k2f car zl;
dlists1:=(car zl . args1) . dlists1>>;
zl:=cdr zl >>;
return dlists1
end;
symbolic procedure exportan dl;
if atom dl then nil
else begin
% extract exp or tan fns from the z-list.
if eq(car dl,'tan) then return t;
nxt: if not eq(car dl,'expt) then return nil;
dl:=cadr dl;
if atom dl then return t;
% Make sure we find nested exponentials?
go to nxt
end;
symbolic procedure findsqrts z;
begin scalar r;
while not null z do <<
if eqcar(car z,'sqrt) then r:=(car z) . r;
z:=cdr z >>;
return r
end;
symbolic procedure trialdiv(x,dl);
begin scalar qlist,q;
while not null dl do
if not null(q:=quotf(x,cdar dl)) then <<
if (caaar dl='tan) and not eqcar(qlist,cdar dl) then
loglist:=('iden . simp cadr caar dl) . loglist;
%tan fiddle!
qlist:=(cdar dl).qlist;
x:=q >>
else dl:=cdr dl;
return qlist.x
end;
endmodule;
module unifac; % Univariate factorization for integration.
% Authors: Mary Ann Moore and Arthur C. Norman.
fluid '(knowndiscrimsign zlist);
global '(!*trint);
exports unifac;
imports cubic,linfac,printdf,quadfac,quadratic,quartic,vp1,
gcd,minusp,prettyprint;
symbolic procedure unifac(pol,var,degree,res);
begin scalar w,q,c;
w:=pol;
if !*trint then superprint w;
%now try looking for linear factors.
trylin: q:=linfac(w);
if null car q then go to nomorelin;
res := ('log . back2df(car q,var)) . res;
w:=cdr q;
go to trylin;
nomorelin:
q:=quadfac(w);
if null car q then go to nomorequad;
res := quadratic(back2df(car q,var),var,res);
w:=cdr q;
go to nomorelin;
nomorequad:
if null w then return res; %all done.
degree:=car w; %degree of what is left.
c:=back2df(w,var);
if degree=3 then res:=cubic(c,var,res)
else if degree=4 then res:=quartic(c,var,res)
else if evenp degree and
pairp (q := halfpower cddr w)
then <<w := (degree/2) . (cadr w . q);
w := unifac(w,var,car w,nil);
res := pluckfactors(w,var,res)>>
else <<
printc "The following has not been split";
printdf c;
res:=('log . c) . res>>;
return res
end;
symbolic procedure halfpower w;
if null w then nil
else if car w=0
then (lambda r;
if r eq 'failed then r else cadr w . r) halfpower cddr w
else 'failed;
symbolic procedure pluckfactors(w,var,res);
begin scalar p,q,knowndiscrimsign;
while w do
<<p := car w;
if car p eq 'atan then nil
else if car p eq 'log
then <<q := doublepower cdr p . q;
%prin2 "q="; %printdf car q;
>>
else interr "Bad form";
w := cdr w>>;
while q do
<<p := car q;
if caaar p=4
then <<knowndiscrimsign := 'negative;
res := quartic(p,var,res);
knowndiscrimsign := nil>>
else if caaar p=2
then res := quadratic(p,var,res)
else res := ('log . p) . res;
q := cdr q>>;
return res
end;
symbolic procedure doublepower r;
if null r then nil
else ((for each j in caar r collect 2*j) . cdar r)
. doublepower cdr r;
symbolic procedure back2df(p,v);
% Undo the effect of uniform.
begin scalar r,n;
n:=car p;
p:=cdr p;
while not minusp n do <<
if not zerop car p then r:=
(vp1(v,n,zlist) .* (car p ./ 1)) .+ r;
p:=cdr p;
n:=n-1 >>;
return reversewoc r
end;
endmodule;
module uniform; % Convert from D.F. to list of coefficients.
% Authors: Mary Ann Moore and Arthur C. Norman.
fluid '(zlist);
exports uniform;
imports exponentof;
symbolic procedure uniform(p,v);
%Convert from d.f. in one variable (v) to a simple list of
%coeffs (with degree consed onto front).
%Fails if coefficients are not all simple integers.
if null p then 0 . (0 . nil)
else begin scalar a,b,c,d;
a:=exponentof(v,lpow p,zlist);
b:=lc p;
if not(denr b=1) then return 'failed;
b:=numr b;
if null b then b:=0
else if not numberp b then return 'failed;
if a=0 then return a . (b . nil); %constant term.
c:=uniform(red p,v);
if c='failed then return 'failed;
d:=car c;
c:=cdr c;
d:=d+1;
while not (a=d) do <<
c:=0 . c;
d:=d+1>>;
return a . (b . c)
end;
endmodule;
module makevars; % Make dummy variables for integration process.
% Authors: Mary Ann Moore and Arthur C. Norman.
fluid '(!*gensymlist!* !*purerisch);
exports getvariables,varsinlist,varsinsq,varsinsf,findzvars,
createindices,mergein;
imports dependsp,union;
% Note that 'i' is already maybe committed for sqrt(-1),
% also 'l' and 'o' are not used as they print badly on certain
% terminals etc and may lead to confusion.
!*gensymlist!* := '(! j ! k ! m ! n ! p ! q ! r ! s ! t ! u ! v ! w ! x
! y ! z);
%mapc(!*gensymlist!*,function remob); %REMOB protection;
symbolic procedure varsinlist(l,vl);
% L is a list of s.q. - find all variables mentioned,
% given thal vl is a list already known about.
begin while not null l do <<
vl:=varsinsf(numr car l,varsinsf(denr car l,vl));
l:=cdr l >>;
return vl
end;
symbolic procedure getvariables sq;
varsinsf(numr sq,varsinsf(denr sq,nil));
symbolic procedure varsinsq(sq,vl);
varsinsf(numr sq,varsinsf(denr sq,vl));
symbolic procedure varsinsf(form,l);
if domainp form then l
else begin
while not domainp form do <<
l:=varsinsf(lc form,union(l,list mvar form));
form:=red form >>;
return l
end;
symbolic procedure findzvars(vl,zl,var,flg);
begin scalar v;
% VL is the crude list of variables found in the original integrand;
% ZL must have merged into it all EXP, LOG etc terms from this.
% If FLG is true then ignore DF as a function.
scan: if null vl then return zl;
v:=car vl; % next variable.
vl:=cdr vl;
% at present items get put onto ZL if they are non-atomic
% and they depend on the main variable. The arguments of
% functions are decomposed by recursive calls to findzvar.
%give up if V has been declared dependent on other things;
if atom v and v neq var and depends(v,var) then
rederr "Can't integrate in the presence of side-relations"
else if not atom v and (not v member zl) and dependsp(v,var)
then if car v='!*sq then zl:=findzvarssq(cadr v,zl,var)
else if car v memq '(times quotient plus minus difference int)
or (((car v) eq 'expt) and fixp caddr v)
then
zl:=findzvars(cdr v,zl,var,flg)
else if flg and car v eq 'df
then <<!*purerisch := t; return zl>> % try and stop it
else zl:=v . findzvars(cdr v,zl,var,flg);
% scan arguments of fn.
%ACH: old code used to look only at CADR if a DF involved.
go to scan
end;
symbolic procedure findzvarssq(sq,zl,var);
findzvarsf(numr sq,findzvarsf(denr sq,zl,var),var);
symbolic procedure findzvarsf(sf,zl,var);
if domainp sf then zl
else findzvarsf(lc sf,
findzvarsf(red sf,
findzvars(list mvar sf,zl,var,nil),
var),
var);
symbolic procedure createindices zl;
% Produces a list of unique indices, each associated with a ;
% different Z-variable;
reversewoc crindex1(zl,!*gensymlist!*);
symbolic procedure crindex1(zl,gl);
begin if null zl then return nil;
if null gl then << gl:=list int!-gensym1 'i; %new symbol needed;
nconc(!*gensymlist!*,gl) >>;
return (car gl) . crindex1(cdr zl,cdr gl) end;
symbolic procedure rmember(a,b);
if null b then nil
else if a=cdar b then car b
else rmember(a,cdr b);
symbolic procedure mergein(dl,ll);
% Adjoin logs of things in dl to existing list ll.
if null dl then ll
else if rmember(car dl,ll) then mergein(cdr dl,ll)
else mergein(cdr dl,('log . car dl) . ll);
endmodule;
module vect; % Vector support routines.
% Authors: Mary Ann Moore and Arthur C. Norman.
% Modified by: James H. Davenport.
exports mkuniquevect,mkvec,mkvecf2q,mkidenm,copyvec,vecsort,swap,
non!-null!-vec,mkvect2;
symbolic procedure mkuniquevect v;
begin scalar u,n;
n:=upbv v;
for i:=0:n do begin
scalar uu;
uu:=getv(v,i);
if not (uu member u)
then u:=uu.u
end;
return mkvec u
end;
symbolic procedure mkvec(l);
begin scalar v,i;
v:=mkvect(isub1 length l);
i:=0;
while l do <<putv(v,i,car l); i:=iadd1 i; l:=cdr l>>;
return v
end;
symbolic procedure mkvecf2q(l);
begin
scalar v,i,ll;
v:=mkvect(isub1 length l);
i:=0;
while l do <<
ll:=car l;
if ll = 0 then ll:=nil;
putv(v,i,!*f2q ll);
i:=iadd1 i;
l:=cdr l >>;
return v
end;
symbolic procedure mkidenm n;
begin
scalar ans,u;
scalar c0,c1;
c0:=nil ./ 1;
c1:= 1 ./ 1;
% constants.
ans:=mkvect(n);
for i:=0 step 1 until n do <<
u:=mkvect n;
for j:=0 step 1 until n do
if i iequal j
then putv(u,j,c1)
else putv(u,j,c0);
putv(ans,i,u) >>;
return ans
end;
symbolic procedure copyvec(v,n);
begin scalar new;
new:=mkvect(n);
for i:=0:n do putv(new,i,getv(v,i));
return new
end;
symbolic procedure vecsort(u,l);
% Sorts vector v of numbers into decreasing order.
% Performs same interchanges of all vectors in the list l.
begin
scalar j,k,n,v,w;
n:=upbv u;% elements 0...n exist.
% algorithm used is a bubble sort.
for i:=1:n do begin
v:=getv(u,i);
k:=i;
loop:
j:=k;
k:=isub1 k;
w:=getv(u,k);
if v<=w
then goto ordered;
putv(u,k,v);
putv(u,j,w);
mapc(l,function (lambda u;swap(u,j,k)));
if k>0
then goto loop;
ordered:
end;
return nil
end;
symbolic procedure swap(u,j,k);
if null u
then nil
else begin
scalar v;
%swaps elements i,j of vector u.
v:=getv(u,j);
putv(u,j,getv(u,k));
putv(u,k,v)
end;
symbolic procedure non!-null!-vec v;
begin
scalar cnt;
cnt := 0;
for i:=0:upbv v do
if getv(v,i)
then cnt:=iadd1 cnt;
return cnt
end;
symbolic procedure mkvect2(n,initial);
begin
scalar u;
u:=mkvect n;
for i:=0:n do
putv(u,i,initial);
return u
end;
endmodule;
end;