module driver; % Driving routines for integration program.
% Author: Mary Ann Moore and Arthur C. Norman.
% Modifications by: John P. Fitch, David Hartley, Francis J. Wright.
fluid '(!*algint
!*backtrace
!*exp
% !*failhard
!*gcd
!*intflag!*
!*keepsqrts
!*limitedfactors
!*mcd
!*noncomp
!*nolnr
!*partialintdf
!*precise
!*purerisch
!*rationalize
!*structure
!*trdint
!*trint
!*uncached
basic!-listofnewsqrts
basic!-listofallsqrts
gaussiani
intvar
kord!*
listofnewsqrts
listofallsqrts
loglist
powlis!*
sqrt!-intvar
sqrt!-places!-alist
subfg!*
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,cflag,intvar,listofnewsqrts,
listofallsqrts,sqrtfn,sqrt!-intvar,sqrt!-places!-alist,
basic!-listofallsqrts,basic!-listofnewsqrts,coefft,
varchange,w,!*precise;
!*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");
then <<varchange := variable . intern gensym();
if !*trint
then printc {"Integration kernel", variable,
"replaced by simple variable", cdr varchange};
variable := cdr varchange>>;
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
<< % added by Alan Barnes
if (cflag:=get(dmode!*, 'cmpxfn)) then onoff('complex, nil);
if (dmod := get(dmode!*,'dname)) then
onoff(dmod,nil)>> where !*msg := nil;
begin scalar dmode!*,!*exp,!*gcd,!*keepsqrts,!*limitedfactors,!*mcd,
!*rationalize,!*structure,!*uncached,kord!*,
ans1,denexp,badbit,nexp,oneterm;
!*keepsqrts := !*limitedfactors := t; % !*sqrt := t;
!*exp := !*gcd := !*mcd := !*structure := !*uncached := t;
dmode!* := nil;
if !*algint
then <<
% The algint code now needs precise off.
% !*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;
if varchange
then <<depend1(car varchange,cdr varchange,t);
expression := int!-subsq(expression,{varchange})>>;
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 := zvars(getvariables x,zv,variable,t);
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);
%% FJW: It is possible for ans here to be just a
%% spurious constant term, in which case we discard it.
if not smemq(variable, ans) then ans := nil ./ 1;
%% This may not be the best place for this fix, but I
%% don't see how it can ever do any harm. [I don't
%% think we need a full depend test here.]
if cdar ans1 neq '(nil . 1)
then ans := addsq(ans,
simpint1(cdar ans1 . variable . w))
>>>>;
end;
ans := multsq(coefft,ans); %Put back coefficient, preserving order.
% 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);
% added by Alan Barnes
if cflag then onoff('complex,t)>> where !*msg := nil;
oldvarstack := varstack!*;
varstack!* := nil;
% ans := errorset!*(list('resimp,mkquote ans),t);
ans := errorset!*(list('int!-resub,mkquote ans,mkquote
varchange),t);
put('int,'simpfn,'simpint);
varstack!* := oldvarstack;
return if errorp ans then error1() else car ans
end;
symbolic procedure int!-resub(x,v);
% {sq,alist} -> sq
% Undo any variable change and resimplify.
if v then <<x := int!-subsq(x,{revpr v}); depend1(car v,cdr v,nil);
resimp x>>
else resimp x;
symbolic procedure int!-subsq(x,v);
% {sq,alist} -> sq
% A version of subsq with the int and df operators unprotected.
% Intended for straightforward change of variable names only.
begin scalar subfuncs,subfg!*;
subfuncs := {remprop('df,'subfunc),remprop('int,'subfunc)};
x := subsq(x,v);
put('df,'subfunc,car subfuncs);
put('int,'subfunc,cadr subfuncs);
return x
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,!*noncomp;
if !*trint then <<
printc "Start of Integration; integrand is ";
printsq integrand >>;
!*noncomp := noncomfp numr integrand
or noncomfp denr integrand;
varlist:=getvariables integrand;
varlist:=varsinlist(xlogs,varlist); %in case more exist in xlogs
if zv then zlist := zv else zlist := zvars(varlist,zlist,var,nil);
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 zvars(x,zv,variable,bool);
% This code attempts to find all possible terms in the target
% integral.
% There used to be problems with nested exponentials or logs,
% but that no longer seems true (10 May 00).
begin scalar oldzlist; integer n;
zv := findzvars(x,list variable,variable,nil);
% The following loop is constrained to five passes to avoid problems
% with differentiation rules such as let {df(f(~x),x) => x*f(x-1)}.
% All integration tests run with just one pass through this loop, so
% five passes is probably overkill.
while oldzlist neq zv and n<5 do <<
oldzlist := zv;
foreach zz in oldzlist do
% zv := findzvars(distexp(pseudodiff(zz,variable)),
% zv,variable,t);
zv := findzvars(pseudodiff(zz,variable),zv,variable,t);
n := n+1>>;
% The following line is based on experiments with the test files.
% At the moment, it's not clear why it's needed, but it is!!
if bool then zv := sort(zv,function ordp);
return zv
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,b,
list('times,nvar,2))), 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))));
if null numr a then return nil; % Power occurred.
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, realdom;
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);
>>;
realdom := not smember('(sqrt -1),integrand);
% print integrand; print realdom;
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 realdom and smember('(sqrt -1),res) then <<
if !*trint then print "Wrong sheet"; 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 symint u % FJW: symbolic integral
end;
mkop 'int!*;
put('int!*,'simpfn,'simpint!*);
symbolic procedure simpint!* u;
begin scalar x;
return if (x := opmtch('int . u)) then simp x
else simpiden('int!* . u)
end;
symbolic procedure remakesf u;
%remakes standard form U, substituting operator INT for INT!*;
if domainp u then u
else addf(multpf(if eqcar(mvar u,'int!*)
then mksp('int . cdr mvar u,ldeg u)
else lpow u,remakesf lc u),
remakesf red u);
symbolic procedure allowedfns u;
if null u then t
else if atom car u 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((((var .^ 1) .* 2) .+ nil) ./ 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.
% FJW: ci,si also defined in specfn(sfint.red), so ...
symbolic((algebraic operator ci,si) where !*msg=nil);
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),
%% FJW: Next 2 rules replaced by ~~ rules below
%% int(e^(~x^2),x) => erf(i*x)*sqrt(pi)/(2i),
%% int(1/e^(~x^2),x) => erf(x) * sqrt(pi)/2,
%% FJW: Missing sqrt(b):
%% int(e^(~b*~x^2),x) => erf(i*x)*sqrt(pi)/(2i*sqrt(b)),
int(e^(~~b*~x^2),x) => erf(i*sqrt(b)*x)*sqrt(pi)/(2i*sqrt(b)),
%% FJW: Rule missing:
int(e^(~x^2/~b),x) => erf(i*x/sqrt(b))*sqrt(pi)*sqrt(b)/(2i),
%% FJW: Missing sqrt(b):
%% int(1/e^(~b*~x^2),x) => erf(x)*sqrt(pi)/(2sqrt(b)),
int(1/e^(~~b*~x^2),x) => erf(sqrt(b)*x)*sqrt(pi)/(2sqrt(b)),
%% FJW: Rule missing:
int(1/e^(~x^2/~b),x) => erf(x/sqrt(b))*sqrt(pi)*sqrt(b)/2,
df(ei(~x),x) => exp(x)/x,
int(e^(~~b*~x)/x,x) => ei(b*x), % FJW
int(e^(~x/~b)/x,x) => ei(x/b),
int(1/(exp(~x*~~b)*x),x) => ei(-x*b), % FJW
int(1/(exp(~x/~b)*x),x) => ei(-x/b),
%% FJW: Next 2 rules replaced by ~~ rules above
%% 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), % FJW
int(sin(~x/~b)/x,x) => si(x/b), % FJW
%% int(sin(~x)/x,x) => si(x), % FJW
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), % FJW
int(cos(~x/~b)/x,x) => ci(x/b), % FJW
%% int(cos(~x)/x,x) => ci(x), % FJW
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, % FJW
int(1/log(~x/~b),x) => ei(log(x/b))*b, % FJW
%% int(1/log(~x),x) => ei(log(x)), % FJW
%% int(1/log(~x+~b),x) => ei(log(x+b)), % FJW
int(1/log(~~a*~x+~b),x) => ei(log(a*x+b))/b, % FJW
int(1/log(~x/~a+~b),x) => ei(log(x/a+b))/b, % FJW
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;
end;