Artifact 4f0e98d75208aa4aa7cb2d3c81a18b4dc5e711b33077538840f7bca0be7c95e6:
- Executable file
r36/src/int.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 154598) [annotate] [blame] [check-ins using] [more...]
module int; % Header for REDUCE integration package. % 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. create!-package('(int contents csolve idepend df2q distrib divide driver intfac ibasics makevars jpatches reform % simplog simpsqrt hacksqrt sqrtf isolve tidysqrt trcase halfangle trialdiv vect dint), % cuberoot d3d4 factr kron lowdeg unifac uniform tdiff '(int trans)); fluid '(!*noextend !*pvar frlis!* gaussiani); global '(gensymcount initl!*); !*pvar:='!_a; gaussiani := !*kk2f '(sqrt -1); gensymcount := 0; initl!* := append('(!*noextend), initl!*); flag('(interr),'transfer); %For the compiler; flag ('(atan dilog ei erf 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. 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 '!!, append(explode u,explode gensymcount))>>; symbolic procedure mknill n; if n=0 then nil else nil . mknill(n-1); symbolic procedure printc u; prin2t u; % This could be an smacro. % Selector written as an smacro. smacro procedure argof u; % Argument of a unary function. cadr 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. endmodule; module contents; % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(content indexlist sqfr varlist zlist); % clogflag 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; for each j in l do res := multf(res,j); % 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; % 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 '(!*trint ccount cmap cmatrix cval loglist neweqn); exports backsubst4cs,createcmap,findpivot,printvecsq, % printspreadc 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(reversip(cvar . reversip 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 Created for log and tan terms:" . 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 idepend; % Routines for considering dependency among variables. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(taylorvariable); exports dependspl,dependsp,involvesq,involvsf; imports 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 smacro procedure taylorfunction u; caar u; 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,w; if null p then return nil ./ 1; d:=denr lc p; w:=red p; while w do <<n:=gcdf(d,denr lc w); % Get denominator of answer ... d:=!*multf(d,quotf(denr lc w,n)); % ... as lcm of denoms in input. w:=red w>>; n := nil; % Place to build numerator of answer. while p do <<w := sqrt2top lc p; n := addf(n,quotf(!*multf(xl2f(lpow p,zlist,indexlist), !*multf(numr w,d)), denr w)); 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, % var2df, 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 VAR2DF(VAR,N,ZLIST); % ((VP1(VAR,N,ZLIST) .* (1 ./ 1)) .+ NIL); 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 needed 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 '(!*trdiv !*trint resid sqrtlist zlist); exports fquotf,testdivdf,dfquotdf; imports df2q,f2df,gcdf,interr,multdf,negdf,plusdf,printdf,printsf, quotf,multsq,invsq,negsq; % Intended for dividing out known factors as produced by the % integration program. Horrible and slow, I expect!! symbolic procedure dfquotdf(a,b); begin scalar resid; if (!*trint or !*trdiv) then << printc "Dividing out a simple factor of "; printdf b >>; a:=dfquotdf1(a,b); if (!*trint or !*trdiv) then << printc "Remaining term to be factorised is "; printdf a >>; if not null resid then begin scalar gres,w; % Make one more check for a zero residue. if null numr df2q resid then return nil; if !*trint or !*trdiv then << printc "Failure in factorisation:"; printdf resid; printc "Which should be zero"; w:=resid; gres:=numr lc w; w:=red w; while not null w do << gres:=gcdf(gres,numr lc w); w:=red w >>; printc "I.e. the following vanishes"; printsf gres>>; interr "Non-exact division due to a log term" end; return a end; symbolic procedure fquotf(a,b); % Input: a and b standard quotients with (a/b) an exact % division with respect to the variables in zlist, % but not necessarily obviously so. the 'non-obvious' problems % will be because of (e.g.) square-root symbols in b % output: standard quotient for (a/b) % (prints message if remainder is not 'clearly' zero. % A must not be zero. begin scalar t1; if null a then interr "a=0 in fquotf"; t1:=quotf(a,b); %try it the easy way if not null t1 then return t1 ./ 1; %ok return df2q dfquotdf(f2df a,f2df b) end; symbolic procedure dfquotdf1(a,b); begin scalar q; if null b then interr "Attempt to divide by zero"; q:=sqrtlist; %remove sqrts from denominator, maybe. while not null q do begin scalar conj; conj:=conjsqrt(b,car q); %conjugate wrt given sqrt if not (b=conj) then << a:=multdf(a,conj); b:=multdf(b,conj) >>; q:=cdr q end; q:=dfquotdf2(a,b); resid:=reversip resid; return q end; symbolic procedure dfquotdf2(a,b); % As above but a and b are distributed forms, as is the result. if null a then nil else begin scalar xd,lcd; xd:=xpdiff(lpow a,lpow b); if xd='failed then << xd:=lt a; a:=red a; resid:=xd .+ resid; return dfquotdf2(a,b) >>; lcd:= !*multsq(lc a,!*invsq lc b); if null numr lcd then return dfquotdf2(red a,b); % Should not be necessary; lcd := xd .* lcd; xd:=plusdf(a,multdf(negdf (lcd .+ nil),b)); if xd and (lpow xd = lpow a % Again, should not be necessary; or xpdiff(lpow xd,lpow b) = 'failed) then <<if !*trint or !*trdiv then <<printc "Problems in dividing:"; printdf xd>>; xd := rootextractdf xd; if !*trint or !*trdiv then printdf xd>>; return lcd .+ dfquotdf2(xd,b) end; symbolic procedure rootextractdf u; if null u then nil else begin scalar v; v := resimp rootextractsq lc u; return if null numr v then rootextractdf red u else (lpow u .* v) .+ rootextractdf red u end; symbolic procedure rootextractsq u; if null numr u then u else rootextractf numr u ./ rootextractf denr u; 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>>>>>>; % D. Dahm suggested an additional call of rootextractf on the % result here. This does cause some expressions to simplify % sooner, but also leads to infinite loops with expressions % like (a*x+b)**p. return if p=0 then addf(c,r) else if null c then r else ((u to p) .* c) .+ r end; % The following hack makes sure that the results of differentiation % gets passed through ROOTEXTRACT % a) This should not be done this way, since the effect is global % b) Should this be done via TIDYSQRT? put('df,'simpfn,'simpdf!*); symbolic procedure simpdf!* u; begin scalar v,v1; v:=simpdf u; v1:=rootextractsq v; if not(v1=v) then return resimp v1 else return v end; symbolic procedure xpdiff(a,b); %Result is list a-b, or 'failed' if a member of this would be negative. if null a then if null b then nil else interr "B too long in xpdiff" else if null b then interr "A too long in xpdiff" else if car b>car a then 'failed else (lambda r; if r='failed then 'failed else (car a-car b) . r) (xpdiff(cdr a,cdr b)); symbolic procedure conjsqrt(b,var); % Subst(var=-var,b). if null b then nil else conjterm(lpow b,lc b,var) .+ conjsqrt(red b,var); symbolic procedure conjterm(xl,coef,var); % Ditto but working on a term. if involvesp(xl,var,zlist) then xl .* negsq coef else xl .* coef; symbolic procedure involvesp(xl,var,zl); % Check if exponent list has non-zero power for variable. if null xl then interr "Var not found in involvesp" else if car zl=var then car xl neq 0 else involvesp(cdr xl,var,cdr zl); endmodule; module driver; % Driving routines for integration program. % Author: Mary Ann Moore and Arthur C. Norman. % Modifications by: John P. Fitch. fluid '(!*algint !*backtrace !*exp !*failhard !*gcd !*intflag!* !*keepsqrts !*limitedfactors !*mcd !*nolnr !*precise !*purerisch !*rationalize !*structure !*trdint !*trint !*uncached basic!-listofnewsqrts basic!-listofallsqrts gaussiani intvar kord!* listofnewsqrts listofallsqrts loglist sqrt!-intvar sqrt!-places!-alist varlist varstack!* xlogs zlist); exports integratesq,simpint,simpint1; imports algebraiccase,algfnpl,findzvars,getvariables,interr,printsq, transcendentalcase,varsinlist,kernp,simpcar,prepsq,mksq,simp, opmtch,formlnr; switch algint,nolnr,trdint,trint; switch hyperbolic; % 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. if atom u or null cdr u or cddr u and (null cdddr u or cddddr u) then rerror(int,1,"Improper number of arguments to INT") else if cddr u then simpdint u % then if getd 'simpdint then simpdint u % else rerror(int,2,"Improper number of arguments to INT") else begin scalar ans,dmod,expression,variable,loglist,oldvarstack, !*intflag!*,!*purerisch,intvar,listofnewsqrts, listofallsqrts,sqrtfn,sqrt!-intvar,sqrt!-places!-alist, basic!-listofallsqrts,basic!-listofnewsqrts,coefft,w; !*intflag!* := t; % Shows we are in integrator. variable := !*a2k cadr u; if not(idp variable or pairp variable and numlistp cdr variable) then typerr(variable,"integration variable"); intvar := variable; % Used in SIMPSQRT and algebraic integrator. w := cddr u; if w then rerror(int,3,"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. % Any changes made to these settings must also be made in wstrass. if dmode!* then onoff(dmod := get(dmode!*,'dname),nil) where !*msg := nil; begin scalar dmode!*,!*exp,!*gcd,!*keepsqrts,!*limitedfactors,!*mcd, !*rationalize,!*structure,!*uncached,kord!*, ans1,denexp,badbit,nexp,oneterm,!*precise; !*keepsqrts := !*limitedfactors := t; % !*sqrt := t; !*exp := !*gcd := !*mcd := !*structure := !*uncached := t; dmode!* := nil; if !*algint then << % The algint code works fine with !*precise on, so use it. !*precise := t; % 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))>>; coefft := (1 ./ 1); % Collect simple coefficients. expression := int!-simp car u; denexp := 1 ./ denr expression; % Get into two bits nexp := numr expression; while not atom nexp and null cdr nexp and not depends(mvar nexp,variable) do <<coefft := multsq(coefft,(((caar nexp) . 1) . nil) ./ 1); nexp := lc nexp>>; ans1 := nil; while nexp do begin % Collect by zvariables scalar x,zv,tmp; if atom nexp then << x:=!*f2q nexp; nexp:=nil >> else << x:=!*t2q car nexp; nexp:=cdr nexp >>; x := multsq(x,denexp); zv := findzvars(getvariables x,list variable,variable,nil); begin scalar oldzlist; while oldzlist neq zv do << oldzlist := zv; foreach zz in oldzlist do zv:=findzvars(distexp(pseudodiff(zz,variable)), zv,variable,t)>>; % The following line was added to make, for example, % int(df(sin(x)/x),x) return the expected result. zv := sort(zv, function ordp) end; tmp := ans1; while tmp do <<if zv=caar tmp then <<rplacd(car tmp,addsq(cdar tmp,x)); tmp := nil; zv := nil>> else tmp := cdr tmp>>; if zv then ans1 := (zv . x) . ans1 end; if length ans1 = 1 then oneterm := t; % Efficiency nexp := ans1; ans := nil ./ 1; badbit:=nil ./ 1; % SQ zero while nexp do % Run down the terms <<u := cdar nexp; if !*trdint then <<princ "Integrate"; printsq u; princ "with Zvars "; print caar nexp>>; ans1 := errorset!*(list('integratesq,mkquote u, mkquote variable,mkquote loglist, mkquote caar nexp), !*backtrace); nexp := cdr nexp; if errorp ans1 then badbit := addsq(badbit,u) else <<ans := addsq(caar ans1, ans); badbit:=addsq(cdar ans1,badbit)>>>>; if !*trdint then <<prin2 "Partial answer="; printsq ans; prin2 "To do="; printsq badbit>>; % We have run down the terms. If there are any bad bits, redo % them. However, since a non-zero badbit implies that % integratesq aborted, the internal variable order may be % confused. So we reset kord!* and reorder expressions in this % case. if badbit neq '(nil . 1) then <<setkorder nil; badbit := reordsq badbit; ans := reordsq ans; coefft := reordsq coefft; if !*trdint then <<princ "Retrying..."; printsq badbit>>; if oneterm and ans = '(nil . 1) then ans1 := nil else ans1 := errorset!*(list('integratesq,mkquote badbit, mkquote variable,mkquote loglist,nil), !*backtrace); if null ans1 or errorp ans1 then ans := addsq(ans,simpint1(badbit . variable . w)) else <<ans := addsq(ans,caar ans1); if cdar ans1 neq '(nil . 1) then ans := addsq(ans, simpint1(cdar ans1 . variable . w)) >>>>; end; ans := multsq(ans, coefft); % Put back coefficient % 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; if !*trdint then << printc "Resimp and all that"; printsq ans >>; % We now need to check that all simplifications have been done % but we have to make sure INT is not resimplified, and that SIMP % does not complain at getting the same argument again. put('int,'simpfn,'simpiden); put('sqrt,'simpfn,sqrtfn); if dmod then onoff(dmod,t) where !*msg := nil; oldvarstack := varstack!*; varstack!* := nil; ans := errorset!*(list('resimp,mkquote ans),t); put('int,'simpfn,'simpint); varstack!* := oldvarstack; return if errorp ans then error1() else car ans end; symbolic procedure numlistp u; % True if u is a list of numbers. null u or numberp car u and numlistp cdr u; % 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 if mvar u = '(sqrt -1) % % then addf(multpf(mksp('i,ldeg 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,zv); begin scalar varlist,x,zlist; if !*trint then << printc "Start of Integration; integrand is "; printsq integrand >>; varlist:=getvariables integrand; varlist:=varsinlist(xlogs,varlist); %in case more exist in xlogs if zv then zlist := zv else <<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 "Determination of the differential field descriptor"; printc "gives the functions:"; print zlist >>; %% Look for rational powers in the descriptor %% If there is make a suitable transformation and do the sub integral %% and return the revised integral x := look_for_substitute(integrand, var, zlist); if x then return x; %% End of rational patch if !*purerisch and not allowedfns zlist then return (nil ./ 1) . integrand; % If it is not suitable for Risch. varlist := setdiff(varlist,zlist); % 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 % **** Treat diffs correctly?? if depends(a,var) then list prepsq simpdf(list(a,var)) else nil else if car a memq '(atan equal log plus quotient sqrt times minus) then begin scalar aa,bb; foreach zz in cdr a do << bb:=pseudodiff(zz,var); aa:= union(bb,aa) >>; return aa end else if car a eq 'expt then if depends(cadr a,var) then if depends(caddr a,var) then prepsq simp list('log,cadr a) . %% a(x)^b(x) cadr a . caddr a . union(pseudodiff(cadr a,var),pseudodiff(caddr a,var)) else cadr a . pseudodiff(cadr a,var) %% a(x)^b else caddr a . pseudodiff(caddr a,var) %% a^b(x) else list prepsq simpdf(list(a,var)); symbolic procedure look_for_substitute(integrand, var, zz); % Search for rational power transformations begin scalar res; if atom zz then return nil else if (res := look_for_rational(integrand, var, zz)) then return res else if (res := look_for_quad(integrand, var, zz)) then return res else if (res := look_for_substitute(integrand, var, car zz)) then return res else return look_for_substitute(integrand, var, cdr zz) end; symbolic procedure look_for_rational(integrand, var, zz); % Look for a form x^(n/m) in the field descriptor, and transform % the integral if it is found. Note that the sqrt form may be used % as well as exponentials. Return nil if no transformation if (car zz = 'sqrt and cadr zz = var) then look_for_rational1(integrand, var, 2) else if (car zz = 'expt) and (cadr zz = var) and (listp caddr zz) and (caaddr zz = 'quotient) and (numberp cadr caddr zz) and (numberp caddr caddr zz) then look_for_rational1(integrand, var, caddr caddr zz) else nil; symbolic procedure look_for_rational1(integrand, var, m); % Actually do the transformation and integral begin scalar newvar, res, ss, mn2m!-1; newvar := gensym(); mn2m!-1 := !*f2q(((newvar .** (m-1)) .* m) .+ nil); %% print ("Integrand was " . integrand); % x => y^m, and dx => m y^(m-1) integrand := multsq(subsq(integrand, list(var . list('expt,newvar,m))), mn2m!-1); if !*trint then << prin2 "Integrand is transformed to "; printsq integrand >>; begin scalar intvar; intvar := newvar; % To circumvent an algint bug. res := integratesq(integrand, newvar, nil, nil); end; ss := list(newvar . list('expt,var, list('quotient, 1, m))); res := subsq(car res, ss) . subsq(quotsq(cdr res, mn2m!-1), ss); if !*trint then << printc "Transforming back..."; printsq car res; prin2 " plus a bad part of "; printsq cdr res >>; return res end; symbolic procedure look_for_quad(integrand, var, zz); % Look for a form sqrt(a+bx+cx^2) in the field descriptor % and transform to the appropriate asin, acosh or asinh. % Return nil if no transformation found if !*algint then nil % as Algint does it better?? else begin if (car zz = 'sqrt and listp cadr zz and caadr zz = 'plus) or (car zz = 'expt and listp cadr zz and caadr zz = 'plus and listp caddr zz and car caddr zz = 'quotient and fixp caddr caddr zz) then << zz := simp cadr zz; if (cdr zz = 1) then << zz := cdr coeff1(prepsq zz, var, nil); if length zz = 2 then return begin % Linear scalar a, b; scalar nvar, res, ss; a := car zz; b := cadr zz; if (depends(a,var) or depends(b,var)) then return nil; nvar := gensym(); if !*trint then << prin2 "Linear shift suggested "; prin2 a; prin2 " "; prin2 b; terpri(); >>; integrand := subsq(integrand, % Make the substitution list(var . list('quotient, list('difference, list('expt,nvar,2),a), b))); integrand := multsq(integrand, % and the dx component simp list('quotient,list('times,nvar,2), b)); % integrand := subsq(integrand, % list(var . list('difference, nvar, a))); % integrand := multsq(integrand, simp b); if !*trint then << prin2 "Integrand is transformed by substitution to "; printsq integrand; prin2 "using substitution "; prin2 var; prin2 " -> "; printsq simp list('quotient, list('difference,list('expt,nvar,2),a), b); >>; res := integratesq(integrand, nvar, nil, nil); ss := list(nvar . list('sqrt,list('plus,list('times,var,b), a))); res := subsq(car res, ss) . subsq(multsq(cdr res, simp list('quotient,2,b)), ss); %% Should one reject if there is a bad bit?? return res; end else if length zz = 3 then return begin % A quadratic scalar a, b, c; a := car zz; b := cadr zz; c:= caddr zz; if (depends(a,var) or depends(b,var) or depends(c,var)) then return nil; a := simp list('difference, a, % Re-centre list('times,b,b, list('quotient,1,list('times,4,c)))); b := simp list('quotient, b, list('times, 2, c)); c := simp c; return if minusf numr c then << if minusf numr a then begin scalar !*hyperbolic; !*hyperbolic := t; return look_for_invhyp(integrand,nil,var,a,b,c) end else look_for_asin(integrand,var,a,b,c)>> else << if minusf numr a then look_for_invhyp(integrand,t,var,a,b,c) else look_for_invhyp(integrand,nil,var,a,b,c) >> end else if length zz = 5 then return begin % A quartic scalar a, b, c, d, e, nn, dd, mm; a := car zz; b := cadr zz; c:= caddr zz; d := cadddr zz; e := car cddddr zz; if not(b = 0) or not(d = 0) then return nil; if (depends(a,var) or depends(c,var)) or depends(e,var) then return nil; nn := numr integrand; dd := denr integrand; if denr(mm :=quotsq(nn ./ 1, !*kk2q var)) = 1 and even_power(numr mm, var) and even_power(dd, var) then << % substitute x -> sqrt(y) return sqrt_substitute(numr mm, dd, var); >>; if denr(mm :=quotsq(dd ./ 1, !*kk2q var)) = 1 and even_power(nn, var) and even_power(numr mm, var) then << % substitute x -> sqrt(y) return sqrt_substitute(nn, multf(dd,!*kk2f var), var); >>; return nil; end; >>>>; return nil; end; symbolic procedure look_for_asin(integrand, var, a, b, c); % Actually do the transformation and integral begin scalar newvar, res, ss, sqmn, onemth, m, n; m := prepsq a; n := prepsq c; b := prepsq b; newvar := gensym(); sqmn := prepsq apply1(get('sqrt, 'simpfn), list list('quotient, list('minus,n), m)); onemth := list('cos, newvar); ss := list('sin, newvar); powlis!* := list(ss, 2, '(nil . t), list('difference,1,list('expt,onemth,2)), nil) . powlis!*; integrand := subs2q multsq(subsq(integrand, list(var . list('difference, list('quotient,ss,sqmn), b))), quotsq(onemth := simp onemth, simp sqmn)); if !*trint then << prin2 "Integrand is transformed by substitution to "; printsq integrand; prin2 "using substitution "; prin2 var; prin2 " -> "; printsq simp list('difference, list('quotient, ss, sqmn), b); >>; res := integratesq(integrand, newvar, nil, nil); powlis!* := cdr powlis!*; ss:= list(newvar . list('asin,list('times,list('plus,var,b),sqmn))); res := subsq(car res, ss) . subsq(quotsq(cdr res, onemth), ss); if !*trint then << printc "Transforming back..."; printsq car res; prin2 " plus a bad part of "; printsq cdr res >>; if (car res = '(nil . 1)) then return nil; return res; end; symbolic procedure look_for_invhyp(integrand, do_acosh, var, a, b, c); % Actually do the transformation and integral; uses acosh/asinh form % depending on second argument begin scalar newvar, res, ss, sqmn, onemth, m, n; m := prepsq a; n := prepsq c; b := prepsq b; newvar := gensym(); if do_acosh then << sqmn := prepsq apply1(get('sqrt, 'simpfn), list list('quotient, n, list('minus, m))); onemth := list('sinh, newvar); ss := list('cosh, newvar) >> else << sqmn:= prepsq apply1(get('sqrt,'simpfn),list list('quotient,n,m)); onemth := list('cosh, newvar); ss := list('sinh, newvar) >>; powlis!* := list(ss, 2, '(nil . t), list((if do_acosh then 'plus else 'difference), list('expt, onemth, 2),1), nil) . powlis!*; % print ("sqmn" . sqmn); print("onemth" . onemth); print ("ss" . ss); % print cdddar powlis!*; integrand := subs2q multsq(subsq(integrand, list(var . list('difference,list('quotient,ss,sqmn),b))), quotsq(onemth := simp onemth, simp sqmn)); if !*trint then << prin2 "Integrand is transformed by substitution to "; printsq integrand; prin2 "using substitution "; prin2 var; prin2 " -> "; printsq simp list('difference, list('quotient, ss, sqmn), b); >>; res := integratesq(integrand, newvar, nil, nil); powlis!* := cdr powlis!*; if !*hyperbolic then << ss := list(if do_acosh then 'acosh else 'asinh, list('times,list('plus,var,b), sqmn)); >> else << ss := list('times,list('plus,var,b), sqmn); ss := if do_acosh then subst(ss,'ss, '(log (plus ss (sqrt (difference (times ss ss) 1))))) else subst(ss,'ss,'(log (plus ss (sqrt (plus (times ss ss) 1))))) >>; ss := list(newvar . ss); res := sqrt2top subsq(car res, ss) . sqrt2top subsq(quotsq(cdr res, onemth), ss); if !*trint then << printc "Transforming back..."; printsq car res; prin2 " plus a bad part of "; printsq cdr res >>; if (car res = '(nil . 1)) then return nil; if smember('(sqrt -1),res) then return nil; % Wrong sheet? return res; end; symbolic procedure simpint1 u; % Varstack* rebound, since FORMLNR use can create recursive % evaluations. (E.g., with int(cos(x)/x**2,x)). begin scalar !*keepsqrts,v,varstack!*; 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; v:=errorset!*(list('simp,mkquote v),!*backtrace); if pairp v then v := car v else v := simp u; !*nolnr := cdr !*nolnr; return v>>; return if (v := opmtch u) then simp v else if !*failhard then rerror(int,4,"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 u 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 then (car u=intvar) or not depends(car u,intvar) else if (caar u = 'expt and not (cadar u = 'e) and not depends(cadar u, intvar) and depends(caddar u, intvar)) then nil else if flagp(caar u,'transcendental) then allowedfns cdr u else nil; symbolic procedure look_for_power(integrand, var); begin scalar nn, dd, mm; nn := numr integrand; dd := denr integrand; if denr(mm :=quotsq(nn ./ 1, !*kk2q var)) = 1 and even_power(numr mm, var) and even_power(dd, var) then << % substitute x -> sqrt(y) return sqrt_substitute(numr mm, dd, var); >>; if denr(mm :=quotsq(dd ./ 1, !*kk2q var)) = 1 and even_power(nn, var) and even_power(numr mm, var) then << % substitute x -> sqrt(y) return sqrt_substitute(nn, numr mm, var); >>; return nil; end; symbolic procedure even_power(xpr, var); if atom xpr then t else if mvar xpr = var then << if evenp pdeg lpow xpr then even_power(lc xpr, var) and even_power(red xpr, var) else nil >> else if eqcar(mvar xpr, 'expt) and cadr mvar xpr = var and evenp caddr mvar xpr then t else if atom mvar xpr then even_power(lc xpr, var) and even_power(red xpr, var) else if even_power(red xpr, var) and even_power(lc xpr, var) then even_prep(mvar xpr, var); symbolic procedure even_prep(xpr,var); if xpr = var then nil else if atom xpr then t else if eqcar(xpr, 'expt) and cadr xpr = var and evenp caddr xpr then t else if even_prep(car xpr, var) then even_prep(cdr xpr, var); symbolic procedure sqrt_substitute(nn, dd, var); begin scalar newvar, integrand, res, ss, !*keepsqrts; newvar := gensym(); integrand := subst(list('sqrt,newvar), var, list('quotient, prepsq (nn ./ dd), 2)); integrand := prepsq simp integrand; integrand := simp integrand; begin scalar intvar; intvar := newvar; % To circumvent an algint bug/oddity res := integratesq(integrand, newvar, nil, nil); end; ss := list(newvar . list('expt, var, 2)); res := subsq(car res, ss) . multsq(list(list(list(var,1),2)) ./ 1, subsq(cdr res, ss)); if !*trint then << printc "Transforming back..."; printsq car res; prin2 " plus a bad part of "; printsq cdr res >>; return res end; % The following rules probably belong in other places. %----------------------------------------------------------------------- algebraic; operator ci,si; % ei. intrules := {e^(~n*acosh(~x)) => (sqrt(x^2-1)+x)^n when numberp n, e^(~n*asinh(~x)) => (sqrt(x^2+1)+x)^n when numberp n, e^(acosh(~x)) => (sqrt(x^2-1)+x), e^(asinh(~x)) => (sqrt(x^2+1)+x), cosh(log(~x)) => (x^2+1)/(2*x), sinh(log(~x)) => (x^2-1)/(2*x), % These next two are rather uncertain. int(log(~x)/(~b-x),x) => dilog(x/b), int(log(~x)/(~b*x-x^2),x) => dilog(x/b)/b + log(x)^2/(2b), int(e^(~x^2),x) => erf(i*x)*sqrt(pi)/(2i), int(1/e^(~x^2),x) => erf(x) * sqrt(pi)/2, int(e^(~b*~x^2),x) => erf(i*x)*sqrt(pi)/(2i*sqrt(b)), int(1/e^(~b*~x^2),x) => erf(x) * sqrt(pi)/(2sqrt(b)), df(ei(~x),x) => exp(x)/x, int(e^(~b*~x)/x,x) => ei(b*x), int(e^(~x/~b)/x,x) => ei(x/b), int(1/(exp(~x*~b)*x),x) => ei(-x*b), int(1/(exp(~x/~b)*x),x) => ei(-x/b), int(e^~x/x,x) => ei(x), int(1/(e^~x*x),x) => ei(-x), int(~a^~x/x,x) => ei(x*log(a)), int(1/((~a^~x)*x),x) => ei(-x*log(a)), df(si(~x),x) => sin(x)/x, int(sin(~b*~x)/x,x) => si(b*x), int(sin(~x)/x,x) => si(x), int(sin(~x)/x^2,x) => -sin(x)/x +ci(x), int(sin(~x)^2/x,x) =>(log(x)-ci(2x))/2, df(ci(~x),x) => cos(x)/x, int(cos(~b*~x)/x,x) => ci(b*x), int(cos(~x)/x,x) => ci(x), int(cos(~x)/x^2,x) => -cos(x)/x -si(x), int(cos(~x)^2/x,x) =>(log(x)+ci(2x)/2), int(1/log(~b*~x),x) => ei(log(b*x))/b, int(1/log(~x),x) => ei(log(x)), int(1/log(~x+~b),x) => ei(log(x+b)), int(1/log(~a*~x+~b),x) => ei(log(a*x+b))/b, int(~x/log(~x),x) => ei(2*log(x)), int(~x^~n/log(x),x) => ei((n+1)*log(x)) when fixp n, int(1/(~x^~n*log(x)),x) => ei((-n+1)*log(x)) when fixp n}; let intrules; endmodule; module intfac; % Interface between integrator and factorizer. % Author: Anthony C. Hearn. % Based on earlier versions by James Davenport, Mary Ann Moore and % Arthur Norman. fluid '(!*intfac !*surds kord!* zlist); % clogflag exports int!-fac; symbolic procedure int!-fac x; % X is a primitive, square-free polynomial, except for monomial % factors. Result is a list of 'factors' wrt zlist. % Throughout most of the integrator we want to add new surds, so we % turn surds on. However, we use *intfac to inhibit use of quadratic % factorizer in the poly/primfac module, since things don't work % properly if this is used. begin scalar !*intfac,!*surds; !*intfac := !*surds := t; return int!-fac!-inner x; end; symbolic procedure int!-fac!-inner x; % X is a primitive, square-free polynomial, except for monomial % factors. Result is a list of 'factors' wrt zlist. begin scalar factors; factors := fctrf x; factors := cdr factors; % Ignore monomial factor. % Make sure x really square-free. factors := for each u in factors collect if cdr u=1 then car u else interr list(x,"not square free"); % It seems we need the logs ordered ahead of atans, hence reverse. return reversip for each u in factors join fac2int u end; symbolic procedure fac2int u; % Returns a list of all the arctangents and logarithms arising from % an attempt to take the one irreducible (but not necessarily the % absolutely irreducible) factor u. begin scalar degrees,x; degrees := for each w in zlist collect (degreef(u,w) . w); if assoc(1,degrees) then return list ('log . (u ./ 1)) % An irreducible polynomial of degree 1 is absolutely irreducible. else if x := assoc(2,degrees) then return int!-quadterm(u,cdr x) else if assoc(0,degrees) then return list('log . (u ./ 1)); % This suggests a surd occurs. Should that be an error? if !*trint then <<printc "*** Polynomial"; printsf u; printc "has not been completely factored">>; return list ('log . (u ./ 1)) end; symbolic procedure int!-quadterm(pol,var); % Add in logs and atans corresponding to splitting the polynomial pol % given it is quadratic wrt var. Does not assume pol is univariate. begin scalar a,b,c,discrim,kord,res,w; kord := setkorder(var . kord!*); % It shouldn't matter if % var occurs twice. c := reorder pol; if ldeg c neq 2 then <<setkorder kord; rerror(int,5,"Invalid polynomial in int-quadterm")>>; a := lc c; c := red c; if not domainp c and mvar c = var and ldeg c = 1 then <<b := lc c; c := red c>>; setkorder kord; discrim := powsubsf addf(multf(b,b),multd(-4,multf(a,c))); if null discrim then interr "discrim is zero in quadterm"; % A quadratic usually implies an atan term. % if not clogflag % then <<w := rootxf(negf discrim,2); % if not(w eq 'failed) then go to atancase>>; w := rootxf(negf discrim,2); if not(w eq 'failed) then go to atancase; w := rootxf(discrim,2); if w eq 'failed then return list ('log . !*f2q pol); % if w eq 'failed then rederr "Integration failure in int-quadterm"; discrim := w; w := multpf(mksp(var,1),a); w := addf(multd(2,w),b); % 2*a*x + b. a := addf(w,discrim); b := addf(w,negf discrim); % Remove monomial multipliers. a := quotf(a,cdr comfac a); b := quotf(b,cdr comfac b); return ('log . !*f2q a) . ('log . !*f2q b) . res; atancase: res := ('log . !*f2q pol) . res; % One part of answer. a := multpf(mksp(var,1),a); a := addf(b,multd(2,a)); a := fquotf(a,w); return ('atan . a) . res end; endmodule; module ibasics; % Some basic support routines for integrator. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(!*backtrace !*gcd !*sqfree !*trint indexlist sqrtflag sqrtlist varlist zlist); exports partialdiff,printdf,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 indx(n); % if n<2 then (list 1) else(n . indx(isub1 n)); symbolic procedure interr mess; <<if !*trint or !*backtrace then <<prin2 "***** INTEGRATION PACKAGE ERROR: "; printc mess>>; error1()>>; 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 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 := setdiff(varlist,zlist); return if sf eq !*sqfree then list list sf else sqfree(sf,zlist) 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 mapply(funct,l); if null l then rerror(int,6,"Empty list to mapply") else if null cdr l then car l else apply2(funct,car l,mapply(funct,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 apply1(f,getv(v,i)) end; endmodule; module makevars; % Make dummy variables for integration process. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(!*gensymlist!* !*purerisch); % exports getvariables,varsinlist,varsinsf,findzvars, % varsinsq % 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 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 % rerror(int,7, % "Can't integrate in the presence of side-relations") zl := union(list v, zl) 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) 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; % printc "Pure set"; 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; reversip 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 cdrmember(a,b); if null b then nil else if a=cdar b then car b else cdrmember(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 cdrmember(car dl,ll) then mergein(cdr dl,ll) else mergein(cdr dl,('log . car dl) . ll); endmodule; module jpatches; % Routines for manipulating sf's with power folding. % Author: James H. Davenport. fluid '(sqrtflag); exports !*addsq,!*multsq,!*invsq,!*multf,!*exptsq,!*exptf; %squashsqrtsq % !*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>>; % In this case (as opposed to ADDSQ), Y and Z may have % developed common factors from SQRT expansion, so a % gcd of Y and Z is needed. x := gcdf(y,z); 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 addf(multpf(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 addf(multpf(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 addf(multpf(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 cdr w neq 0 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 rerror(int,8,"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)) % Cannot be *MULTF. 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 (if null x then squashsqrt red sf else lpow sf .* x .+ squashsqrt red sf) where x = squashsqrt lc sf; %remd 'simpx1; % The following definition requires frlis!* declared global. %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 intersection(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); % Raises 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); % Raises 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 reform; % Reformulate expressions using C-constant substitution. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(!*trint cmap cval loglist ulist); exports logstosq,substinulist; imports prepsq,mksp,nth,multsq,addsq,domainp,invsq,plusdf; symbolic procedure substinulist ulst; % Substitutes for the C-constants in the values of the U's given in % ULST. Result is a D.F. if null ulst then nil else begin scalar temp,lcu; lcu:=lc ulst; temp:=evaluateuconst numr lcu; if null numr temp then temp:=nil else temp:=((lpow ulst) .* !*multsq(temp,!*invsq(denr lcu ./ 1))) .+ nil; return plusdf(temp,substinulist red ulst) 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 simpsqrt; % Simplify square roots. % Authors: Mary Ann Moore and Arthur C. Norman. % Heavily modified by J.H. Davenport for algebraic functions. fluid '(!*galois !*pvar !*tra !*trint basic!-listofallsqrts gaussiani basic!-listofnewsqrts intvar kord!* knowntobeindep listofallsqrts listofnewsqrts sqrtflag sqrtlist sqrt!-places!-alist varlist zlist); % 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 doesn't have an AVALUE then it is itself. if car sq eq 'times then return mapply(function multsq, for each j in cdr sq collect simpsqrti j); 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); w:=errorset!*(list('afactor,mkquote v,mkquote !*pvar),t); 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 := mksqrt u; listofnewsqrts:=w . listofnewsqrts; return !*kk2f w; notnew: w:=car w; v:=stt(w,!*pvar); if car v neq 1 then errach list("Error in mknewsqrt: ",v); 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 errach list("Division failure in mknewsqrt",u); return w end; symbolic procedure actualsimpsqrt(sf); if sf iequal -1 then gaussiani else actualsqrtinner(sf,listofnewsqrts); symbolic procedure actualsqrtinner(sf,l); if sf =1 then 1 else if null l or domainp sf or ldeg sf=1 % Patch by A.C. Norman to prevent recursion errors. then actualsimpsqrt2 sf else begin scalar z; if numberp sf and (z := list('sqrt,sf)) member l then return !*kk2f z; z := argof car l; if z member l then z := !*kk2f car l else z := !*q2f simp z; if z = -1 then return actualsqrtinner(sf,cdr l); z:=quotf(sf,z); if null z then return actualsqrtinner(sf,cdr l) else 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 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 intvar and involvesf(sf,intvar) then intvar else findatom mvar sf); % intvar is the best thing to do square-free % decompositions with respect to, but anything % else will do if intvar 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; 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:= for each u in things collect sublis(u,thisplace) . u; pos := for each u in answers collect (thisplace . car u) . cdr u; % pos is sqrt(f) -> sqrt(innersubst f) neg := for each u in answers collect {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 sqrtf; % Square root of standard forms. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(!*noextend zlist); exports nrootn,domainp,minusf; % minusdfp,sqrtdf 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 isolve; % Routines for solving the final reduction equation. % Author: Mary Ann Moore and Arthur C. Norman. % Modifications by: John P. Fitch. fluid '(!*trint badpart ccount cmap cmatrix cval indexlist lhs!* lorder orderofelim power!-list!* pt rhs!* sillieslist tanlist ulist zlist); global '(!*number!* !*statistics); 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). begin top: if null lhs!* then return ulist else begin scalar u,lpowlhs; lpowlhs:=lpow lhs!*; begin scalar ll,m1,chge; ll:=maxorder(power!-list!*,zlist,0); m1:=lorder; while m1 do << if car ll < car m1 then << chge:=t; rplaca(m1,car ll) >>; ll:=cdr ll; m1:=cdr m1 >>; if !*trint and chge then << princ "Maximum order for undetermined coefficients is reduced to "; printc lorder >> end; u:=pickupu(rhs!*,lpow lhs!*,t); if null u then << if !*trint then << printc "***** Equation for a constant to be solved:"; 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!*); u := t >> 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 << printc "A coefficient of numerator has been determined"; 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 and u then << printc "Terms remaining are:"; printdf lhs!*; printc " ">>; end; go to top 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 constant "; prin2 subby; prin2 " by "; printsq simp xp >>; sillieslist:=subby . sillieslist; return subdf(express,xp,subby) end; symbolic procedure checku(ulst,u); % Checks that U is not already in ULST - ie. that this u-coeff % has not already been given a value. ulst and (car u = caar ulst or checku(cdr ulst,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. pu and ((car pu+caar pr)<0 or 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,mesgiven,needsquash; 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 << if not mesgiven then << %% Changed by JPff prin2t "*** Introduce new constants for coefficients"; mesgiven := t >>; 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 terms: "; prin2 car xx; princ "*("; printsq cdr xx; printc ")">>; if (not testord(car xx,lorder)) then << if !*trint then 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 << needsquash := t; 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); if needsquash then lhs!*:=squashconstants(lhs!*) >>; 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 tansfrom(rhs,z,indexlist,n); % We notice that in all bad cases we have (j-num)tan**j...; % Extract the num to get list of all maxima; if null z then nil else begin scalar zz,r, rr, ans; r:=rhs; zz := car z; ans := 0; if not(atom zz) and car zz = 'tan then while r do << rr:=caar r; % The list of powers; for i:=1:n do rr:=cdr rr; if fixp caar rr then ans := max(ans,tanextract(car indexlist,prepsq cdar r)); r:=cdr r; >>; return cons(ans,tansfrom(rhs, cdr z,cdr indexlist,n+1)) end; symbolic procedure tanextract(var, exp); % Find the value of the variable which makes the expression vanish. % The coefficients must be linear. begin scalar ans, c0, c1; ans := cdr coeff1(exp,var,nil); if length ans = 2 and not(car ans = 0) then << c0 := car ans; c1 := cadr ans; if eqcar(c0,'!*sq) then c0 := cadr c0 else c0 := c0 ./ 1; if eqcar(c1,'!*sq) then c1 := cadr c1 else c1 := c1 ./ 1; ans := multsq(c0, invsq c1); if atom ans then return 0; if (cdr ans = 1) and fixp (car ans) then return -(car ans); return 0 >>; return 0; 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 constant 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(minpowers,z,n); % Find a limit on the order of terms, this is ad hoc; if null z then nil else if eqcar(car z,'sqrt) then cons(1,maxorder(cdr minpowers,cdr z,n+1)) else if (atom car z) or (caar z neq 'tan) then cons(maxfrom(lhs!*,n)+1,maxorder(cdr minpowers,cdr z,n+1)) else cons(max(car minpowers, maxfrom(lhs!*,n)), maxorder(cdr minpowers,cdr z,n+1)); symbolic procedure maxfrom(l,n); maxfrom1(l,n+1,0); symbolic procedure maxfrom1(l,n,v); % Largest order in the nth variable. if null l then v else <<v := max(nth(caar l,n),v); maxfrom1(cdr l,n,v)>>; 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 := for each j in cmatrix collect addtovector j 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 reversip 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 ccount=0 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 "This equation adds no new information"; return nil>>; % equation was 0=0. noway: badpart:=tokill . badpart; % non-integrable term. if !*trint then <<printc "Inconsistency in equations for constants,"; printc " so non integrable">>; 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 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 n1,dd; n1:=tidysqrtf numr q; if null n1 then return nil ./ 1; %answer is zero. dd:=tidysqrtf denr q; return multsq(n1,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 cdr x neq 0 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 !*failhard !*nowarnings !*purerisch !*reverse !*trint badpart ccount cmap cmatrix content cval denbad denominator!* indexlist lhs!* loglist lorder orderofelim power!-list!* rhs!* sillieslist sqfr sqrtflag sqrtlist tanlist svar varlist xlogs zlist); % !*reverse: flag to re-order zlist. % !*nowarnings: flag to lose messages. global '(!*number!* !*ratintspecial !*seplogs !*spsize!* !*statistics 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; % 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,power!-list!*, sqrtlist,tanlist,loglist,dflogs,eprim,dfun,unintegrand, sqrtflag,badpart,rhs!*,lhs!*,gcdq,cmap,cval,orderofelim,cmatrix; scalar 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); % subs2q? 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)); % multf? if !*trint then << printc "After unnormalization the integrand is "; printsq unintegrand >> end; divlist := findtrialdivs zlist; % Also puts some things on loglist sometimes. sqrtlist := findsqrts zlist; divlist := trialdiv(denr unintegrand,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:"; % prettyprint sqfr>>; if !*reverse then zlist := reverse zlist; % Alter order function. indexlist := createindices zlist; % if !*trint then << printc "...indices are:"; % prettyprint indexlist>>; dfu:=dfnumr(svar,car divlist); % if !*trint then << terpri(); % printc "************ Derivative of u is:"; % printsq dfu>>; loglist := append(loglist,factorlistlist prim); %%% nconc? 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 <<if !*trint then printsq checkdffail(rhs!*,svar); interr "Simplification fails on above expression">>; if !*trint then << printc "Distributed Form of Numerator is:"; printdf rhs!*>>; lhs!* := f2df dflogs; % if checkdffail(lhs!*,svar) then interr "Simplification failure"; if !*trint then << printc "Distributed Form of integrand is:"; printdf lhs!*; terpri()>>; cval := mkvect(ccount); for i := 0:ccount do putv(cval,i,nil ./ 1); power!-list!* := tansfrom(rhs!*,zlist,indexlist,0); lorder:=maxorder(power!-list!*,zlist,0); originalorder := for each x in lorder collect x; % Must copy as it is overwritten. if !*trint then << printc "Maximum order for variables 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 ccount neq 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; result:= !*multsq(dfun,!*invsq(denominator!* ./ 1)); result:= !*multsq(result,!*invsq(jhd!-content ./ 1)); 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 << %% prettyprint result; terpri(); printc "*****************************************************"; printc "************ THE INTEGRAL IS : **********************"; printc "*****************************************************"; terpri(); printsq result; terpri()>>; if badpart then begin scalar n,oorder; if !*trint then printc "plus a part which has not been integrated"; lhs!*:=badpart; lorder:=maxorder(power!-list!*,zlist,0); oorder:=originalorder; n:=length lorder; while lorder do << if car lorder > car originalorder then wrongway:=t; if car lorder=car originalorder then n:= n-1; lorder:=cdr lorder; originalorder:=cdr originalorder >>; %% if n=0 then wrongway:=t; % Nothing changed if !*trint and wrongway then printc "Went wrong way"; dfun:=df2q badpart; %% if !*trint %% then <<printsq dfun; printc "Denbad = "; printsf denbad>>; dfun:= !*multsq(dfun,invsq(denbad ./ 1)); badpart := dfun; %%% ***** if wrongway then << if !*trint then printc "Resetting...."; result:=nil ./ 1; dfun := integrand; badpart:=dfun >>; if rootcheckp(unintegrand,svar) then return simpint1(integrand . svar.nil) . (nil ./ 1) else if !*purerisch or allowedfns zlist then << badpart := dfun; dfun := nil ./ 1 >> % JPff else << !*purerisch:=t; if !*trint then <<printc " Applying transformations ..."; printsq dfun>>; denbad:=transform(dfun,svar); if denbad=dfun then << dfun:=nil ./ 1; badpart:=denbad >> else <<denbad:=errorset!*(list('integratesq,mkquote denbad, mkquote svar,mkquote xlogs, nil), !*backtrace); %%% JPF fix for distributive version if not atom denbad then << denbad:=car denbad; % As from an errorset dfun:=untan car denbad; if (dfun neq '(nil . 1)) then badpart:=untan cdr denbad; % There is still a chance that we went the wrong way. % Decide that it is better if there is no bad part if car badpart and not(badpart=denbad) then << wrongway:=nil; lhs!*:=f2df car badpart; lorder:=maxorder(power!-list!*,zlist,0); n:=length lorder; while lorder do << if car lorder > car oorder then wrongway:=t; if car lorder=car oorder then n:= n-1; lorder:=cdr lorder; oorder:=cdr oorder >>; if wrongway or (n=0) then << if !*trint then printc "Still backwards"; dfun := nil ./ 1; badpart := integrand>>>>>> else <<badpart := dfun; dfun := nil ./ 1 >>>>>>; if !*failhard then rerror(int,9,"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) end else badpart:=nil ./ 1; return (sqrt2top result . badpart) % JPff end; symbolic procedure checkdffail(u,v); % Sometimes simplification fails and this gives the integrator the % idea that something is a constant when it is not. Check for this. if null u then nil else if depends(lc u,v) then lc u else 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 "; prin2 i; prin2t ":"; terpri()>>; 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. % Do this with regard to cases like sin(x)*tan(x/2), so attempt to % limit times we use half angles. begin scalar zl,tnarg,substlist; zl := zlist; while car(tnarg := tan!-function!-in zl) and halfangle!-confusion(zlist,cadar tnarg) do <<zl := cdr tnarg; tnarg := car tnarg; if eqcar(tnarg,'tan) then substlist := (gensym() . tnarg) . substlist else substlist := (gensym() . list('quotient,1,('tan . cdr tnarg))) . substlist; u := subst(caar substlist,tnarg,u)>>; return if substlist then simp sublis(substlist,prepsq halfangle(u,x)) % simp prepsq was added so that 1/(e**x*cos(1/e**x)**2) % for example returns a closed-form result. else simp prepsq halfangle(u,x) end; symbolic procedure tan!-function!-in zz; % Look at zlist for tangents or cotangents. <<while zz and not eqcar(car zz,'tan) and not eqcar(car zz,'cot) do zz := cdr zz; if null zz then nil . nil else zz>>; symbolic procedure halfangle!-confusion(zz,tnarg); % Is there a function in the zlist with twice the tangent argument? <<while zz and (atom car zz or not(tnarg = prepsq simp list('quotient,cadar zz,2))) do zz := cdr zz; zz>>; 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); if cadr k eq 'pi then return !*k2q k; % Don't consider tan(pi/2). 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) = 'cot then quotqq(!*subtrq(!*int2qm(1), exptsq(kt,2)),multsq(!*int2qm(2),kt)) else if car(k) = 'sec then quotqq(addsq(!*int2qm(1), exptsq(kt,2)), !*subtrq(!*int2qm(1),exptsq(kt,2))) else if car(k) = 'csc then quotqq(addsq(!*int2qm(1),exptsq(kt,2)), !*subtrq(!*int2qm(1),exptsq(kt,2))) else if car(k) = 'sinh then quotqq(!*subtrq(!*p2q mksp('expt.('e. cdr k),2), !*int2qm(1)), multsq(!*int2qm(2), !*p2q mksp('expt . ('e . cdr(k)),1))) else if car(k) = 'cosh then quotqq(addsq(!*p2q mksp('expt.('e. cdr k),2), !*int2qm(1)), multsq(!*int2qm(2), !*p2q mksp('expt . ('e . cdr(k)),1))) else if car(k) = 'tanh then quotqq(!*subtrq(!*p2q mksp('expt.('e. cdr k),2), !*int2qm(1)), addsq(!*p2q mksp ('expt.('e.cdr(k)),2), !*int2qm(1))) else if car(k) = 'coth then quotqq(addsq(!*p2q mksp('expt.('e.cdr(k)),2), !*int2qm(1)), !*subtrq(!*p2q mksp('expt.('e . cdr k),2),!*int2qm(1))) else if car(k) = 'acot then !*p2q mksp(list('atan, list('quotient, 1, cadr k)),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. % We turn off gcd to avoid unnecessary gcd calculations, as suggested % by Rainer Schoepf. begin scalar !*gcd,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 begin scalar yy; %% printc "Recursive tan"; printc cadr y; yy := prepsq untan simp cadr y . nil; %% princ "==> "; printc yy; if evenp z then <<z := z/2; return simp list('quotient, list('plus, list('minus, list('cos, 'times . (2 . yy))), 1),list('plus, list('cos, 'times . (2 . yy)), 1))>> else if z=1 then return simp list('quotient, list('plus, list('minus, list('cos, 'times . (2 . yy))), 1),list('sin, 'times . (2 . yy))) else <<z := (z - 1)/2; w := simp list('quotient, list('plus, list('minus, list('cos, 'times . (2 . yy))), 1),list('sin, 'times . (2 . yy))); return simp list('quotient, list('plus, list('minus, list('cos, 'times . (2 . yy))), 1),list('plus, list('cos, 'times . (2 . yy)), 1)) >> end 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 '(!*trint loglist tanlist); 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 reversip 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 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; module dint; % Definite integration support. % Author: Anthony C. Hearn. fluid '(!*precise); symbolic procedure simpdint u; begin scalar low,upp,fn,var,x,y; if length u neq 4 then rerror(int,2,"Improper number of arguments to INT"); fn := car u; var := cadr u; low := caddr u; upp := cadddr u; low := reval low; upp := reval upp; if low = upp then return nil ./ 1 else if null getd 'new_defint then nil else if upp = 'infinity then if low = 0 then if not smemql('(infinity unknown), x := defint!* {fn,var}) then return simp!* x else nil else if low = '(minus infinity) then return mkinfint(fn,var) else if freeof(var,low) then if not smemql('(infinity unknown), x := defint!* {fn,var}) and not smemql('(infinity unknown), y := indefint!* {fn,var,low}) then return simp!* {'difference,x,y} else nil else nil else if upp = '(minus infinity) or low = 'infinity then return negsq simpdint {fn,var,upp,low} else if low = '(minus infinity) then return simpdint{prepsq simp{'sub,{'equal,var,{'minus,var}},fn}, var,{'minus,upp},'infinity} else if low = 0 then if freeof(var,upp) and not smemql('(infinity unknown), x := indefint!* {fn,var,upp}) then return simp!* x else nil else if freeof(var,upp) and freeof(var,low) and not smemq('(infinity unknown), x := indefint!* {fn,var,upp}) and not smemql('(infinity unknown), y := indefint!* {fn,var,low}) then return simp!* {'difference,x,y}; return mkdint(fn,var,low,upp) end; symbolic procedure defint!* u; (if errorp x then 'unknown else car x) where x = errorset2 {'new_defint,mkquote u}; symbolic procedure indefint!* u; (if errorp x or eqcar(car x,'indefint2) then 'unknown else car x) where x = errorset2 {'new_indefint,mkquote u}; symbolic procedure mkdint(fn,var,low,upp); % This could be used as an entry point to other dint procedures. % Should we handle infinity, - infinity differently? begin scalar x,!*precise; if getd 'defint0 and not((x := defint0 {fn,var,low,upp}) eq 'failed) then return simp x else if not smemq('infinity,low) and not smemq('infinity,upp) then <<x := prepsq!* simpint {fn,var}; if not eqcar(x,'int) then return simp!* {'difference, subeval {{'equal,var,upp},x}, subeval {{'equal,var,low},x}}>>; return mksq({'int,fn,var,low,upp},1) end; symbolic procedure mkinfint(fn,var); begin scalar x,y; if getd 'defint0 and not((x := defint0 {fn,var,'(minus infinity),'infinity}) eq 'failed) then return simp x; x := simpdint {fn,var,0,'infinity}; y := simpdint {fn,var,'(minus infinity),0}; if kernp x and eqcar(mvar numr x,'int) and kernp y and eqcar(mvar numr y,'int) then return mkdint(fn,var,'(minus infinity),'infinity) else return addsq(x,y) end; endmodule; end;