module symint; % Improved simplification of symbolic integrals
% Author: Francis J. Wright
% An extension of simpint1 in module driver (by Mary Ann Moore, Arthur
% C. Norman and John P. Fitch) to provide better simplification of
% integrals of symbolic derivatives and integrals. (Originally
% motivated by the needs of the CRACK package.)
% Change Log:
% 7/1/98: Partial integration for integrals of integrals
% 10/1/98: Extended partial integration for integrals of integrals
% 11/1/98: Commutation of integrals
% 21/2/98: df(y,x,2) etc. handling corrected
fluid '(!*failhard !*IntDfFound);
switch CommuteInt; % off by default (for now)
deflist('((CommuteInt ((t (rmsubs))))), 'simpfg);
% If the switch CommuteInt is turned on then the top-level integration
% in a symbolic integral is commuted into the integrand to try to
% simplify it, and if that fails and the result is a symbolic multiple
% integral then it is left in canonical nesting order (as is already
% done automatically for multiple derivatives). However, an
% integrable nested derivative is integrated regardless of this
% switch.
switch PartialInt, PartialIntDf, PartialIntInt; % off by default
deflist('((PartialInt ((t (on '(PartialIntDf PartialIntInt)))
(nil (off '(PartialIntDf PartialIntInt)))))
(PartialIntDf ((t (rmsubs))))
(PartialIntInt ((t (rmsubs))))), 'simpfg);
% If the switch PartialIntDf is turned on then integration by parts is
% performed if the result simplifies in the sense that it integrates a
% symbolic derivative and does not introduce new symbolic derivatives.
% However, because the initial integral contains an unevaluated
% derivative then the result must still contain an unevaluated
% integral.
% If the switch PartialIntInt is turned on then integration by parts
% is performed if the result simplifies in the sense that it removes a
% symbolic integral from the integrand and does not introduce new
% symbolic integrals. However, because the initial integral contains
% an unevaluated integral then the result must still contain an
% unevaluated integral.
% The switch PartialInt is just a convenience to turn both the above
% switches on or off together.
switch XPartialInt, XPartialIntDf, XPartialIntInt; % off by default
deflist('((XPartialInt ((t (on '(XPartialIntDf XPartialIntInt)))
(nil (off '(XPartialIntDf XPartialIntInt)))))
(XPartialIntDf ((t (rmsubs))))
(XPartialIntInt ((t (rmsubs))))), 'simpfg);
% These switches control extended partial integration of integrals of
% the form int( int(u(x,z),z) * v(x), x ), which is experimental,
% somewhat heuristic and may be slow.
symbolic procedure symint u;
% u has the form (int y x).
% At this point linearity has been applied.
begin scalar v, y, x;
y := cadr u; x := caddr u;
% Check for a directly integrable derivative:
if (v := NestedIntDf(y,x,nil)) then return mksq(v,1);
if !*failhard then rerror(int,4,"FAILHARD switch set");
if (!*PartialIntDf or !*PartialIntInt) and
% Integrate by parts if the result simplifies:
% DO WE NEED TO CALL SIMPINT1 RECURSIVELY ON THE RESULT?
(v := PartialInt(y,x)) then return mksq(v,1);
if (!*XPartialIntDf or !*XPartialIntInt) and
% EXPERIMENTAL! Try extended partial integration:
(v := XPartialInt(y,x)) then return mksq(v,1);
return mksq(u,1)
end;
%% symbolic procedure NestedIntDf(y, x);
%% %% int( ... df(f,A,x,B) ..., x) -> ... df(f,A,B) ...
%% %% Find a df(f,A,x,B) among possibly nested int's and df's within
%% %% the integrand y in int(y,x), and return the whole structure y
%% %% but with the derivative integrated; otherwise return nil.
%% %% [A,B are arbitrary sequences of kernels.]
%% not atom y and
%% begin scalar car_y, nested;
%% return
%% if (car_y := car y) eq 'df and memq(x, cddr y) then
%% %% int( df(f, A, x, B), x ) -> df(f, A, B)
%% 'df . cadr y . delete(x, cddr y)
%% %% use delete for portability!
%% %% deleq is defined in CSL, delq in PSL -- oops!
%% else if memq(car_y, '(df int)) and
%% (nested := NestedIntDf(cadr y, x)) then
%% %% int( df(int(df(f, A, x, B), c), C), x ) ->
%% %% df(int(df(f, A, B), c), C)
%% %% int( int(df(f, A, x, B), c), x ) ->
%% %% int(df(f, A, B), c)
%% car_y . nested . cddr y
%% end;
symbolic procedure NestedIntDf(y, x, !*recursive);
%% In order to simplify a symbolic integral int(y,x), commute the
%% integral through integrals and derivatives in the integrand to
%% try to find an integrable integrand. Return the result if
%% successful; otherwise return nil. [A,B are arbitrary sequences
%% of kernels or "kernels followed by integers".] If the integral
%% does not simplify, optionally commute multiple integrals into
%% canonical nesting order, as is done in the standard
%% differentiator code for multiple derivatives. !*recursive is
%% nil in the top-level call, t in recursive calls. [NB: The
%% top-level call of this procedure makes redundant the let rule in
%% the standard integrator code to integrate derivatives.]
not atom y and
begin scalar fn, nested;
return
if (fn := car y) eq 'df then % integrating a derivative
if (nested := IntDf(y, x)) then nested
%% int( ... df(f, A, x, B) ... , x ) -> df(f, A, B)
else if (nested := NestedIntDf(cadr y, x, t)) then
%% recursing into the integrand
fn . nested . cddr y
else nil
else if !*failhard then nil % give up!
else if fn eq 'int then % integrating an integral
if eq(x, caddr y) then
%% int( ... int(f, x) ... , x ) -> stop
nil
else if (nested := NestedIntDf(cadr y, x, t)) then
%% recursing into the integrand
fn . nested . cddr y
else if !*CommuteInt and ordp(x, caddr y) then
%% Commute integrals into canonical nesting order:
%% int( ... int(f, b) ... , a ) ->
%% int( ... int(f, a) ... , b )
%% Successive calls of the integrator by the simplifier
%% to integrate nested integrals causes this code to
%% sort the integrands into canonical order.
{'int, {'int, cadr y, x}, caddr y}
else nil
else if !*recursive and !*CommuteInt and
%% y is not an integral or a derivative -- try to
%% integrate it unless at top level:
not eqcar(nested := reval {'int,y,x}, 'int) then nested
end;
symbolic procedure IntDf(y, x);
% y = df(f, u, nu, v, nv, ...) where nu, nv, ... optional
% if x = u, v, ... then return int(y, x)
begin scalar !*IntDfFound;
x := IntDfVars(cddr y, x);
if !*IntDfFound then return
if x then 'df . cadr y . x else cadr y
end;
symbolic procedure IntDfVars(y, x);
if y then
if car y eq x then
begin scalar n;
!*IntDfFound := t;
return
if (y := cdr y) and fixp(n := car y) then <<
y := cdr y;
if n > 2 then y := (n-1) . y;
x . y
>> else y
end
else car y . IntDfVars(cdr y, x);
symbolic procedure PartialInt(y, x);
%% Integrate by parts if the resulting integral simplifies;
%% otherwise return nil. Split integrand into a derivative or
%% integral and a second factor and call the appropriate procedure.
%% Try all possible allowed partial integrations in turn.
not atom y and
begin scalar denlist, faclist, facs, df_or_int, result;
% Process any quotient:
if car y eq 'quotient then <<
denlist := cddr y;
% y := numerator:
if atom(y := cadr y) then return % no derivative or integral
>>;
% y := list of factors:
if car y eq 'times then y := cdr y
else if denlist or !*PartialIntInt then y := y . nil
% Can do double integral int(int(u(x),x),x) as a special case
else return;
faclist := y;
% Loop through all integrable derivatives or differentiable
% integrals among the factors:
continue:
while faclist and ( atom(df_or_int := car faclist) or
not (memq(car df_or_int, '(df int)) and
memq(x, cddr df_or_int)) ) do
faclist := cdr faclist;
% Finally, break the loop if there is no integrable derivative
% or differentiable integral:
if null faclist then return;
facs := delete(df_or_int, y); % list of factors
facs := if null facs then 1
else if cdr facs then 'times . facs else car facs;
if denlist then facs := 'quotient . facs . denlist;
if car df_or_int eq 'df then
(if !*PartialIntDf and
(result := PartialIntDf(facs, df_or_int, x))
then return result)
else
(if !*PartialIntInt and
(result := PartialIntInt(df_or_int, facs, x))
then return result);
% Continue the loop through the factors in faclist:
faclist := cdr faclist;
goto continue
end;
symbolic procedure PartialIntDf(u, df_v, x);
%% int(u(x)*df(v(x),x), x) -> u(x)*v(x) - int(df(u(x),x)*v(x), x)
%% Integrate by parts if the resulting integral simplifies [to
%% avoid infinite loops], which means that df(u(x),x) may not
%% contain any unevaluated derivatives; otherwise return nil.
begin scalar v;
v := IntDf(df_v, x);
% Check that df(u(x),x) simplifies:
if smemq('df, df_v := reval {'df,u,x}) then return;
return reval {'difference,
{'times,u,v}, {'int, {'times, df_v, v}, x}}
end;
symbolic procedure PartialIntInt(int_u, v, x);
%% int(int(u(x),x) * v(x), x) ->
%% int(u(x),x) * int(v(x),x) - int( u(x) * int(v(x),x), x )
%% Integrate by parts if the resulting integral simplifies [to
%% avoid infinite loops], which means that int(v(x),x) may not
%% remain an unevaluated integral; otherwise return nil.
begin scalar u;
u := cadr int_u; % kernel being integrated
% Check that int(v(x),x) simplifies:
if eqcar(v := reval {'int,v,x}, 'int) then return;
return reval {'difference,
{'times,int_u,v}, {'int, {'times,u,v}, x}}
end;
symbolic procedure XPartialInt(y, x);
%% Extended partial integration. This code is somewhat heuristic
%% and may be slow. The problem is to try to simplify
%% int( int(u(x,z),z) * v(x), x ).
%% Integrate by parts if the resulting integral simplifies;
%% otherwise return nil. Split integrand into an integral NOT wrt
%% x and a second factor and call the appropriate procedure. Try
%% all possible allowed partial integrations in turn.
not atom y and
begin scalar denlist, faclist, facs, int, result;
% Process any quotient:
if car y eq 'quotient then <<
denlist := cddr y;
% y := numerator:
if atom(y := cadr y) then return % no derivative or integral
>>;
% y := list of factors:
if car y eq 'times then y := cdr y
else if denlist then y := y . nil
else return;
faclist := y;
% Loop through all integrals among the factors:
continue:
while faclist and ( not eqcar(int := car faclist, 'int) or
eq(x, caddr int) ) do
faclist := cdr faclist;
% Finally, break the loop if there is no appropriate integral:
if null faclist then return;
facs := delete(int, y); % list of factors
facs := if null facs then 1
else if cdr facs then 'times . facs else car facs;
if denlist then facs := 'quotient . facs . denlist;
if facs = 1 then return; % ?????
if (result :=
(!*XPartialIntDf and XPartialIntDf(int, facs, x)) or
(!*XPartialIntInt and XPartialIntInt(int, facs, x)))
then return result;
% Continue the loop through the factors in faclist:
faclist := cdr faclist;
goto continue
end;
symbolic procedure XPartialIntDf(int_u, v, x);
%% int(int(u(x,z),z)*v(x), x) ->
%% int(int(u(x,z),x),z) * v(x) -
%% int( int(int(u(x,z),x),z) * df(v(x),x), x )
%% provided df(v(x),z) and int(u(x,z),x) simplify.
begin scalar df_v, z;
% Check that df(v(x),x) simplifies:
if smemq('df, df_v := reval {'df, v, x}) then return;
z := caddr int_u;
int_u := reval {'int, cadr int_u, x}; % int(u(x,z),x)
% Check that int(u(x,z),x) simplifies:
if eqcar(int_u, 'int) then return;
int_u := reval {'int, int_u, z}; % int(int(u(x,z),x),z)
return reval {'difference,
{'times,int_u,v}, {'int, {'times, int_u, df_v}, x}}
end;
symbolic procedure XPartialIntInt(int_u, v, x);
%% int(int(u(x,z),z) * v(x), x) ->
%% int(u(x,z),z) * int(v(x),x) -
%% int( int(df(u(x,z),x),z) * int(v(x),x), x )
%% provided int(v(x),x) and int(df(u(x,z),x),z) simplify.
%% If x = z then this reduces to PartialIntInt.
begin scalar u;
u := cadr int_u; % kernel being integrated
% Check that int(v(x),x) simplifies:
if eqcar(v := reval {'int, v, x}, 'int) then return;
% Check that df(u(x),x) simplifies:
if smemq('df, u := reval {'df, u, x}) then return;
% Check that int(df(u(x,z),x),z) simplifies:
if eqcar(u := reval {'int, u, caddr int_u}, 'int) then return;
return reval {'difference,
{'times,int_u,v}, {'int, {'times,u,v}, x}}
end;
endmodule;
end;