Artifact 88ff4b895761900bc16ae0a6e47f36d213704d7139f73a84fb5e0dc804523094:
- Executable file
r37/packages/defint/defintx.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: 24842) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/defint/defintx.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: 24842) [annotate] [blame] [check-ins using]
module defintx; % Code for definite integration using contour methods. % Author: Stanley L. Kameny <stan_kameny@rand.org>. load_package solve,misc; fluid '(!*allpoly rdon!* !*norationalgi); switch allpoly; global '(domainlist!* poles!*); algebraic << logcomplex := { log(~x + i) => log(sqrt(x*x+1))+i*atan2(1/sqrt(x*x+1),x/sqrt(x*x+1)) when repart(x)=x, log(~x - i) => log(sqrt(x*x+1))-i*atan2(1/sqrt(x*x+1),x/sqrt(x*x+1)) when repart(x)=x, log(~x + i*~y) => log(sqrt(x*x+y*y))+i*atan2(y/sqrt(x*x+y*y),x/sqrt(x*x+y*y)) when repart(x)=x and repart(y)=y, log(~x - i*~y) => log(sqrt(x*x+y*y))-i*atan2(y/sqrt(x*x+y*y),x/sqrt(x*x+y*y)) when repart(x)=x and repart(y)=y, log(~x/~y) => log(x) - log(y) when repart(y)=y, log(sqrt ~x) => (log x)/2, log(-1) => i*pi, log (-i) => -i*pi/2, log i => i*pi/2, log(-~x) => i*pi+log x when repart(x)=x and numberp x and x>0, log(-i*~x) => -i*pi/2 + log x when repart(x)=x and numberp x and x>0, log(i*~x) => i*pi/2 + log x when repart(x)=x and numberp x and x>0 }$ atan2eval := { atan2(sqrt 3/2,1/2) => pi/3, atan2(-sqrt 3/2,1/2) => -pi/3, atan2(sqrt 3/2,-1/2) => 2*pi/3, atan2(-sqrt 3/2,-1/2) => -2*pi/3, atan2(3/(2*sqrt 3),1/2) => pi/3, % these shouldn't be needed atan2(-3/(2*sqrt 3),1/2) => -pi/3, % these shouldn't be needed atan2(3/(2*sqrt 3),-1/2) => 2*pi/3, % these shouldn't be needed atan2(-3/(2*sqrt 3),-1/2) => -2*pi/3, % these shouldn't be needed atan2(1/2,sqrt 3/2) => pi/6, atan2(-1/2,sqrt 3/2) => -pi/6, atan2(1/2,-sqrt 3/2) => 5*pi/6, atan2(-1/2,-sqrt 3/2) => -5*pi/6, atan2(1/2,3/(2*sqrt 3)) => pi/6, % these shouldn't be needed atan2(-1/2,3/(2*sqrt 3)) => -pi/6, % these shouldn't be needed atan2(1/2,-3/(2*sqrt 3)) => 5*pi/6, % these shouldn't be needed atan2(-1/2,-3*(2*sqrt 3)) => -5*pi/6, % these shouldn't be needed atan2(sqrt 2/2,sqrt 2/2) => pi/4, atan2(-sqrt 2/2,sqrt 2/2) => -pi/4, atan2(sqrt 2/2,-sqrt 2/2) => 3*pi/4, atan2(-sqrt 2/2,-sqrt 2/2) => -3*pi/4, atan2(0,-1) => pi, atan2(0,1) => 0, atan2(1,0) => pi/2, atan2(-1,0) => -pi/2 }$ ipower := {i^~n => cos(n*pi/2) + i*sin(n*pi/2), (-i)^~n => cos(n*pi/2) - i*sin(n*pi/2)}$ >> $ algebraic let ipower,atan2eval; %algebraic let logcomplex,atan2eval; fluid '(!*diffsoln zplist!! poles!# !*msg !*rounded !*complex zlist); switch diffsoln; load_package int; % put('defint,'psopfn,'defint0); symbolic procedure defint0 u; begin scalar rdon!*,!*msg,c,!*noneglogs,fac,!*norationalgi, !*combinelogs,!*rationalize; if not getd 'solvesq then load_package solve; if length u neq 4 then rederr "defint called with wrong number of args"; c := !*complex; off complex; % since complex on won't work here! % on complex; % this causes trouble here, so it was moved into % defint11s after splitfactors has operated! !*noneglogs := t; algebraic (let logcomplex); %,atan2eval); fac := !*factor; on factor; !*norationalgi := t; u := errorset2 {'defint1,mkquote u}; !*norationalgi := nil; if errorp u then <<u := 'failed; go to ret>> else u := car u; off factor; if !*rounded then % if approximate answer, eliminate infinitessimal real or % imaginary part. (<<off complex; if domainp numr u and denr u = 1 then (if evallessp({'times,{'abs,prepsq im},eps}, {'abs,prepsq rl}) then u := rl else if evallessp({'times,{'abs,prepsq rl},eps}, {'abs,prepsq im}) then u := addsq(u,negsq rl)); u := mk!*sq u; if rdon!* then off rounded;off complex; go to ret2>> where rl=repartsq u,im=impartsq u,eps=10.0^(2-precision 0)); !*rationalize := t; u := aeval prepsq u; on complex; u := simp!* u; % u := evalletsub2({'(logcomplexs), % {'simp!*,{'prepsq,mkquote u}}},nil); % if errorp u then error(99,list("error during log simp")) % else u := car u; ret: if fac then on factor; algebraic (clearrules logcomplex); %,atan2eval); if u neq 'failed then u := prepsq u; off complex; on combinelogs; if u neq 'failed then u := aeval u; ret2: if c then on complex; return u end; symbolic procedure defint1 u; defint11s(car u,cadr u,caddr u,cadddr u); symbolic procedure defint11s(exp,var,llim,ulim); % removes from integrand any factors independent of var, and passes % the dependent factors to defint11. Based on FACTOR being on. <<% off complex; % not needed here since complex is off already. exp := splitfactors(simp!* aeval exp,var); on complex; % at this point it is safe to turn complex on. multsq(simp!* car exp, defint11(cadr exp,var,llim,ulim,t))>>; symbolic procedure fxinfinity x; if x eq 'infinity then 'inf else if x = '(minus infinity) then 'minf else x; symbolic procedure defint11(exp,var,llim,ulim,dtst); if evalequal(llim := fxinfinity llim, ulim := fxinfinity ulim) or evalequal(exp,0) then nil ./ 1 else begin scalar !*norationalgi,r,p,q,poles,rlrts,cmprts,q1; scalar m,n; if ulim = 'minf or llim = 'inf then return defint11({'minus,exp},var,ulim,llim,dtst); exp := simp!* exp; % Now the lower limit must be 'minf or a finite value, % and the upper limit must be 'inf or a finite value. There are % four cases: % Upper limit is 'inf. Convert lower limit to zero if necessary, % and use methods for infinite integrals. if ulim = 'inf then <<if not(llim member '(0 minf)) then <<exp := subsq(exp,{var . {'plus,var,llim}}); llim := 0>>; go to c0>>; % lower limit is 'minf. Convert this case to upper limit 'inf. if llim = 'minf then <<off complex; exp := reval prepsq subsq(exp,{var . {'minus,var}}); llim := reval {'minus,ulim}; on complex; return defint11(exp,var,llim,'inf,dtst)>>; % Both limits are finite, so check for indef integral and % substitute values if it exists; else check for special forms with % finite limits, try substitutions, or abort. r := simpint {prepsq exp,var}; if eqcar(prepsq r,'int) then go to c1; p := errorset2 list('subsq, mkquote r, mkquote {var . ulim}); q := errorset2 list('subsq, mkquote r, mkquote {var . llim}); if errorp(p) or errorp (q) then << p:= simplimit list('limit!- ,mk!*sq(r),var,ulim); q:= simplimit list('limit!+ ,mk!*sq(r),var,llim); >> else <<p := car p; q := car q>>; return q1 := addsq(p,negsq q); c1: rederr "special forms for finite limits not implemented"; c0: r := exp; p := numr r; q := denr r; % if not polynomp(q,var) then % rederr "only polynomial denominator implemented"; m := degreeof(p,var); n := degreeof(q,var); if smemql('(failed infinity),m) or smemql('(failed infinity),n) then return error(99, 'failed); % Note that degreeof may return a fraction or a general complex % quantity. if not evalgreaterp(prepsq addsq(repartsq n,negsq repartsq m),1) then go to div; % this is the point at which special cases can be tested. if (q1 := specformtestint(q,p,var,llim,ulim)) then return q1; % beyond here, only rational functions are handled. if not (m := sq2int m) or not (n := sq2int n) then <<write "this irrational function case not handled"; terpri(); error(99,'failed)>>; if n - m < 2 then go to div; if dtst and !*diffsoln then if (q1 := diffsol(q,p,m,n,var,llim,ulim)) then return q1; off factor; !*norationalgi := nil; poles := getpoles(q,var,llim); rlrts := append(car poles,cadr poles); cmprts := caddr poles; !*norationalgi := t; q1 := difff(q,var); q := q ./ 1; p := p ./ 1; return if llim = 0 then defint2(p,q,q1,var,rlrts,cmprts) else defint3(p,q,q1,var,rlrts,cmprts); div: % write "integral diverges"; terpri(); error(99,'failed) end; symbolic procedure zpsubsq x; subsq(x,for each v in zplist!! collect (v . 0)); symbolic procedure degreeof(p,var); % p is a standard form. % Note that limit returns "failed" as a structure, not an id. % Also, the limit code has problems with bessels at the present time. % if smemq('besseli,p) then !*k2q 'failed else if smemql ('(besselj besselk bessely besseli),p) then !*k2q 'failed else (if null car de then de else <<if d then onoff(d := get(d,'dname),nil); p := simp!* limit(list('quotient,list('times,var, prepsq de),prepf p), var,'infinity); if d then onoff(d,t); p>>) where d=dmode!*,de=difff(p,var); symbolic procedure genminusp x; if domainp x then !:minusp x else !:minusp topeval prepf x; symbolic procedure sq2int x; (if null numr impartsq x and denr y=1 then if null z then 0 else if numberp z then z else nil) where z=numr y where y=repartsq x; symbolic procedure topeval u; <<if not r then on rounded; if not c then on complex; u := numr simp!* aeval u; if not r then off rounded;if not c then off complex; u>> where r=!*rounded,c=!*complex,!*msg=nil; symbolic procedure firstatom x; if atom x then x else firstatom car x; symbolic procedure valueof u; (if firstatom x neq 'root_of then x) where x=caar u; symbolic procedure rdsolvesq u; solvesq(subf(numr simp!* cadr x,list((caddr x) . caadr u)), caadr u,caddr u) where x=caaaar caar u; symbolic procedure defint2(p,q,q1,var,rlrts,cmprts); % Does the actual computation of integral with limits 0, inf. % Pertinent poles and their orders have been found previously. begin scalar int; p := simp!* aeval{'times,{'log,{'minus,var}},prepsq p}; int := nil ./ 1; for each r in append(rlrts,cmprts) do int := addsq(int,residuum(p,q,q1,var,prepsq car r,cdr r)); return negsq int end; symbolic procedure defint3(p,q,q1,var,rlrts,cmprts); % Does the actual computation of integral with limits minf, inf. % Pertinent poles and their orders have been found previously. begin scalar int,int2; int := int2 := nil ./ 1; for each r in cmprts do int := addsq(int,residuum(p,q,q1,var,prepsq car r,cdr r)); int := addsq(int,int); for each r in rlrts do int2 := addsq(int2,residuum(p,q,q1,var,prepsq car r,cdr r)); int := addsq(int,int2); return multsq(simp!* '(times pi i),int) end; symbolic procedure diffsqn(sq,var,n); <<if n>0 then for j := 1:n do sq := diffsq(sq,var); sq>>; symbolic procedure polypwrp(exp,var); begin scalar pol,fl; integer s,pwr; if eqcar(exp,'expt) then <<pol := cadr exp; if (pwr := caddr exp) <2 then return nil; if atom pol then if var eq pol then s := 1 else return nil else if not eqcar(pol,'plus) then return nil else for each p in cdr pol do s := max(s,termvpwr(p,var)); return if s = 0 then nil else {pol,s,pwr}>> else if eqcar(exp,'times) then <<exp := for each p in cdr exp collect polypwrp(p,var); for each p in exp do <<if null p then fl := t; if not fl then pwr := gcdn(pwr,caddr p)>>; if fl then return nil; s := (for each p in exp sum cadr p * caddr p)/pwr; pol := 'times . for each p in exp collect {'expt,car p,caddr p/pwr}; return {pol,s,pwr}>> end; symbolic procedure termvpwr(p,var); if freeof(p,var) then 0 else if atom p then 1 else if eqcar(p,'expt) and cadr p = var and numberp caddr p then caddr p else if eqcar(p,'times) then for each q in cdr p sum termvpwr(q,var) else 0; symbolic procedure diffsol(q,p,mm,nn,var,llim,ulim); % p is numerator q is denom mm is deg p nn is deg q (q := polypwrp(prepf q,var)) and begin scalar n,s,m,r,zplist!!; n := mm; s := cadr q; m := caddr q; % if s, the power of the base polynomial, > 4 then the base % polynomial won't be solved, and this approach won't work! % However, for s > 2, the approach is impractical, because the % roots of the zp!! polynomial are too complicated, so in the % following, s is tested s > 2. if s > 2 or m*s neq nn or nn - n <= 2 then return nil; r := (n+2)/s; if r*s < n+2 then r := r+1; if m = r then return nil; q := {'plus,car q,'zp!!}; zplist!! := '(zp!!); q := numr simp!*{'expt,q,r}; nn :=(-1)^(m - r)*factorial(r - 1) ./ factorial(m - 1); p := defint11(prepsq(p ./ q),var,llim,ulim,nil); p := zpsubsq diffsqn(p,'zp!!,m - r); return multsq(nn,p) end; symbolic procedure residuum(p,q,q1,var,pole,m); if m=1 then subsq(quotsq(p,q1),list(var . pole)) else begin integer n; q1 := nil; for each r in poles!* do <<n := cdr r; r := prepsq car r; if not evalequal(pole,r) then q1 := {'expt,{'difference,var,r},n} . q1>>; n := ((lc numr simp!* prepsq q) where !*factor=nil); q1 := 'times . (n . q1); return if ((lt numr simp!* prepsq q = lt numr simp!*{'times,{'expt,{'difference,var,pole},m},q1}) where !*factor=nil) then <<q := quotsq(p,simp!* q1); q := diffsqn(q,var,m - 1); q := subsq(q,{var . pole}); q := if null numr q then q else quotsq(q,factorial(m -1) ./ 1)>> else q1 := simp!* (p := limit( prepsq quotsq(diffsqn(multsq(quotsq(p,q), simp!* {'expt,{'difference,var,pole},m}),var,m - 1), factorial(m - 1) ./ 1),var,pole)) end; symbolic procedure splitfactors(u,var); % returns a list of two factors: % independent of var and dependent on var. begin scalar n,d,ni,nd,di,dd,fli,fld; n := prepf numr u; if n=0 then return {0,0}; d := prepf denr u; ni := nd := di := dd := 1; if simptermp n then <<if freeof(n,var) then ni := n else nd := n; go to d>>; for each x in cdr n do if freeof(x,var) then ni := if ni = 1 then list x else <<fli := t; x . ni>> else nd := if nd = 1 then list x else <<fld := t; x . nd>>; ni := fleshoutt(ni,fli); nd := fleshoutt(nd,fld); fli := fld := nil; d: if simptermp d then <<if freeof(d,var) then di := d else dd := d; go to ret>>; for each x in cdr d do if freeof(x,var) then di := if di = 1 then list x else <<fli := t; x . di>> else dd := if dd = 1 then list x else <<fld := t; x . dd>>; di := fleshoutt(di,fli); dd := fleshoutt(dd,fld); ret: return {fleshout(ni,di),fleshout(nd,dd)} end; symbolic procedure simptermp x; atom x or ((y = 'minus and simptermp cadr x or y neq 'times) where y=car x); symbolic procedure fleshout(n,d); if d = 1 then n else {'quotient,n,d}; symbolic procedure fleshoutt(n,d); if n = 1 then n else if d then 'times . n else car n; symbolic procedure specformtestint(den,num,var,llim,ulim); % This tests for defint(x^(p-1)/(a*x^n+b)^m,x,0,inf) with % m,n,p positive integers, p/n not integer and m>(p/n) and either % a*b>0 or {a*b<0,m=1}. % Since splitfactors has removed all factors which do not depend upon % var, both num and den are either 1 or products of terms which % depend upon var. begin scalar a,b,d,m,n,p,q1,q,k,z,ff; den := prepf den; num := prepf num; if not(llim=0) and ulim='inf then go to t2; % This is the test for defint(y**(q-1)/(a*y**n +b)**m,y,0,inf); % which is converted to defint(x**(p-1)/(x+z)**m,x,0,inf); % the next test is assumed to be accessed at label t2. if num = 1 then q1 := 0 else if (q1 := varpwrtst(num,var))=0 then go to t2; if atom den then go to t2 else if not eqcar(den,'times) then %only (a*y**n+b)**m term in den. if (k := aynbmtst(den,var)) then go to sep4 else go to t2 else if length den neq 3 then go to t2; % the denominator is the product of 2 terms, one of which must be % y**q and the other an aynbm form like the previous case. d := cdr den; if not((k := aynbmtst(cadr d,var)) and eqcar(q := varpwrtst(car d,var),'quotient) or (k := aynbmtst(car d,var)) and eqcar(q := varpwrtst(cadr d,var),'quotient)) then go to t2; sep4: n := caddr k; if not numberp n then go to t3; q := prepsq simp!* {'plus,1,q1,{'minus,q}}; p := prepsq simp!* {'quotient,q,n}; m := cadddr k; if not numberp m or m<1 then go to t3; a := car k; b := cadr k; z := prepsq simp!* {'quotient,b,a}; if numr impartsq simp!* z then go to t2; ff := prepsq simp!* aeval {'quotient,1,{'times,n,{'expt,a,m}}}; % there are two different cases: % z > 0 and m >repart p >0 m >= 1 % z < 0 and m=1 (Cauchy principal value) if evalgreaterp(z,0) then if not (evalgreaterp((k := prepsq repartsq simp!* p),0) and evalgreaterp(m,k)) then go to t3 else <<k := prepsq simp!* aeval {'times,{'expt,-1,m+1},'pi,{'expt,z,{'difference,p,m}}}; n := 1; for c := 1:(m-1) do n := prepsq simp!* aeval {'times,n,{'difference,p,c}}; q := prepsq simp!* aeval {'times,{'factorial,m-1},{'sin,{'times,p,'pi}}}; return simp!* aeval {'quotient,{'times,k,n,ff},q}>>; if m neq 1 then go to t3; write "Cauchy principal value"; terpri(); k := prepsq simp!* aeval {'minus,{'expt,{'quotient,-1,z},{'difference,1,p}}}; q := prepsq simp!* aeval {'times,ff,{'quotient,'pi,{'times,a,n}},{'cot,{'times,'pi,p}}}; return simp!* aeval {'times,k,q}; t3: return nil; % most (if not all) of these are divergence cases. t2: return specformtestint2(den,num,var,llim,ulim) end; symbolic procedure aynbmtst(exp,var); % test for form (a*y**n+b)**m (or degenerate forms of this) and % extract parameters a,n,b,m. b qnd m are required to be present. % car exp = 'expt or else m=1. begin scalar a,b,m,n; if not eqcar(exp,'expt) then <<m := 1; goto a2>>; m := caddr exp; exp := cadr exp; a2: if not eqcar(exp,'plus) or length exp neq 3 then return nil; b := caddr exp; if eqcar(cadr exp,'times) then <<exp := cdadr exp; if length exp neq 2 or not( numberp (a := car exp) and (n := varpwrtst(cadr exp,var)) neq 0 or numberp (a := cadr exp) and (n := varpwrtst(car exp,var)) neq 0) then return nil>> else <<a := 1; if (n := varpwrtst(cadr exp,var)) = 0 then return nil>>; return {a,b,n,m} end; fluid '(!*fullpoly); switch fullpoly; symbolic procedure getpoles(q,var,llim); begin scalar poles,rt,m,rlrt,cmprt,rtv,rtvz,cpv,prlrts,nrlrts,rlrts, cmprts,!*multiplicities,!*fullroots,!*norationalgi; off factor; !*norationalgi := poles!* := nil; !*multiplicities := t; if !*fullpoly then !*fullroots := t; % if !*allpoly = 'all then % <<on rounded; rdon!* := t; write "test mode"; terpri()>>; poles := solvesq(simp!* prepf q,var,1); !*norationalgi := t; lp: if null poles then go to int; rt := car poles; poles := cdr poles; m := caddr rt; rlrt := cmprt := nil; if (rtv := valueof rt) then <<poles!* := (rtv . m) . poles!*; rtvz := zpsubsq rtv; rt := car impartsq rtvz; if null rt or (rt := topevalsetsq prepf rt) and evalequal(0,prepsq rt) then rlrt := rtv else cmprt := rtv; if llim = 0 then <<if rlrt then <<if null car rtvz then go to div else if not genminusp car rtvz then <<if m > 1 then go to div else cpv := t; prlrts := (rlrt . m) . prlrts>> else nrlrts := (rlrt . m) . nrlrts>> else cmprts := (cmprt . m) . cmprts; go to lp>> else <<if rlrt then <<if m > 1 then go to div else cpv := t; rlrts := (rlrt . m) . rlrts>> else if not genminusp car impartsq rtvz then cmprts := (cmprt . m) . cmprts>>; go to lp>>; una: if !*rounded then rederr "unable to find poles approximately"; if not !*allpoly then <<write "Denominator degree > 4. Approx integral requires ON ALLPOLY"; terpri(); error(99,"failed")>> else <<write "approximate integral only"; terpri()>>; on rounded; rdon!* := t; if valueof car(rt := rdsolvesq rt) then <<poles := append(rt,poles); go to lp>>; go to una; div: % write "integral diverges"; terpri(); error(99,'failed); int: if cpv then <<write "Cauchy principal value"; terpri()>>; return if llim=0 then {prlrts,nrlrts,cmprts} else {rlrts,nil,cmprts} end; symbolic procedure specformtestint2(den,num,var,llim,ulim); % This checks for defint(x^k*R(x),x,0 inf) where {k != 0,-1<k<1} and % limit+(x^(k+1)*R(x),x,0)=limit(x^(k+1)*R(x),x,inf)=0 where R is % a rational function with no poles of order >1 on positive real axis. begin scalar d,k,k1,m,p,p1,q,q1,poles,prpoles,s1,s2; if not(llim=0) and ulim='inf then go to t2; p1 := polanalyz(num,var); k1 := polanalyz(den,var); if not (p1 or k1) then go to t2; k := prepsq simp!* aeval {'difference,p1,k1}; if numberp k or not(evalgreaterp(k,-1) and evalgreaterp(1,k)) then go to t2;%<== this was t3 but caused problem! if (d := dmode!*) then onoff(d := get(d,'dname),nil); m := prepsq simp!* aeval {'quotient,{'times,var,num},den}; if numr simp!* limit!+(m,var,0) or numr simp!* limit(m,var,'infinity) then go to t3; if d then onoff(d,t); % all tests met, except for finding poles of den. % move irrational factor to numerator, changing the sign of var. p := simp!* aeval {'times,num, {'expt,var,{'times,-1,p1}},{'expt,{'minus,var},k}}; % note that p in general has a non-trivial denominator. % now remove irrational factor from denominator, leaving polynomial. q := simp!* aeval {'times,den,{'expt,var,{'times,-1,k1}}}; q1 := diffsq(q,var); poles := getpoles(numr q,var,llim); prpoles := car poles; poles := append(cadr poles,caddr poles); s1 := s2 := nil ./ 1; k1 := {'times,'pi,{'plus,k,1}}; if poles then <<for each r in poles do s1 := addsq(s1,residuum(p,q,q1,var,prepsq car r,cdr r)); s1 := {'quotient,{'times,'pi,prepsq s1},{'sin,k1}}>> else s1 := 0; if prpoles then <<for each r in prpoles do s2 := addsq(s2,residuum(p,q,q1,var,prepsq car r,cdr r)); s2 := {'times,'pi,prepsq s2,{'cot,k1}}>> else s2 := 0; return simp!* aeval {'difference,s1,s2}; t2: return nil; % replace by call to next test. t3: % write "integral diverges"; terpri(); error(99,'failed) end; symbolic procedure polanalyz(exp,var); % test for fractional power of var in exp. if poltstp exp then ((if eqcar( exp := varpwrtst(if eqcar(ex2,'times) then cadr ex2 else ex2,var), 'quotient) then exp else 0) where ex2=if eqcar(exp,'minus) then cadr exp else exp); symbolic procedure poltstp exp; atom exp and exp or car exp member domainlist!* or car exp member '(plus times quotient minus expt sqrt) and begin scalar fg; for each c in cdr exp do if not poltstp c then fg := t; return null fg end; symbolic procedure evalmax(a,b); if numberp a and numberp b then max(a,b) else if evalgreaterp(a,b) then a else b; symbolic procedure evalplus(a,b); if numberp a and numberp b then a+b else prepsq simp!* aeval {'plus,a,b}; symbolic procedure varpwrtst(exp,var); if atom exp then if exp = var then 1 else 0 else if car exp eq 'minus then varpwrtst(cadr exp,var) else if car exp member '(plus,difference) then (<<for each c in cdr exp do q := evalmax(q,varpwrtst(c,var)); q>> where q=0) else if eqcar(exp,'expt) then prepsq simp!* aeval{'times,varpwrtst(cadr exp,var),caddr exp} else if eqcar(exp,'sqrt) then prepsq simp!* aeval{'times,varpwrtst(cadr exp,var),{'quotient,1,2}} else if eqcar(exp,'times) then (<<for each c in cdr exp do q := evalplus(q,varpwrtst(c,var)); q>> where q=0) else 0; endmodule; end;