reduce-historical/r37/packages/defint/defintx.red

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;