\documentclass{article}
\usepackage{fullpage,amsmath}
\title{DRAFT:
Partitioning of algebraic subexpressions:
in Computer Algebra Systems.
Don't use matching so much.
An integration example too.
}
\author{Richard Fateman}
\begin{document}
\maketitle
\begin{abstract}
A popular technique to direct transformations of algebraic expressions
in a computer algebra system such as Macsyma/Maxima or Mathematica is
to write a pattern and use a system routine to match it, followed by
some rule-replacement. This works nicely for many applications, but
is inefficient and sometimes ineffective, especially if the pattern
must match an expression headed by n-ary commutative operator with
identity such as ``+'' or ``*''. That is, we wish to avoid
unification with identity, and the backup it often entails: compounding
the issue is that writing correct programs can be difficult.
We demonstrate an alternative. We especially show how this can work
nicely for writing a symbolic integration program; a table-driven
half-page of code that does many standard problems.
\end{abstract}
\section{Introduction}
It happens that programmers wish to partition sums or products into
pieces that match various predicates. For example, consider an
integration program that tries to work up a case for solving an
expression like $ax^3\cos(x)$. This is not such a problem if we know
in advance the symbolic value of $x$, say $t$. Then $\cos(t)$ can be
divided out and we can search for $t^n$, and look for its coefficient,
$a$ which is presumably some expression that is free of $t$. Maxima's
pattern matcher is pretty good at this. But less so in trying
something like this: $f(u)u'$, that is, where one item is a derivative
of something embedded in the other, and $f$ could be any of a number
of functions.
These are tough problems for ``pattern matching'', but such problems
don't have to be done that way.
\section {Partitioner}
We describe a program that works on sums, products, or perhaps other
expressions with a different head operator. It is given a predicate
{\tt p}, a program that returns true or false given an
expression. Roughly speaking it separates the given expression into
pieces that pass and fail the predicate. To make the program useful,
we set up 4 other arguments, described below, and then we have a
fairly general partitioner.
Since we will collect those pieces that pass the predicate and those
that don't, how shall we collect them? Let $y$ (for ``yes'') be the
collection that passes the predicate. Initialize $y$ to a parameter
{\tt init} which might be 1 for products, 0 for sums. For each term
$v$ for which {\tt p($v$)} is true, we combine $v$ with $y$ by a {\tt
combiner} for example, we usually find it convenient to combine using
``*'' in the case of products, ``+'' for sums. Similarly, $n$ (for
``no'').
We are almost done: what do we do with the collected $y$ and $n$? We
could return them as a pair, or try something more devious that would
work with other matching facilities. Here's what we programmed. Let
the user of the partition program provide 2 more items. The first is
an {\tt action} which could just be a dummy function name, or even
``[`` to return a list of [$y$,$n$]. The second is a symbolic name
{\tt res}, whose value will be bound to the result. This final item is
something of a kludge, but makes it easier for certain restricted
constructions in pattern matching to get a handle on the result.
\begin{verbatim}
/* Our partitioning tool. Note that the expression being partitioned
is the last argument.
This makes matchdeclare simpler */
partition_expression(operator,pred,init,combiner,action,res,E):=
block([yes:init,no:init],
if not atom(E) and inpart(E,0)=operator then
(map(lambda([r], if apply(pred,[r])=false
then no:apply(combiner,[r,no])
else yes:apply(combiner,[r,yes])),
inargs(E) ),
res :: apply(action, [yes,no]), /* result stored as requested */
true)),
inargs(z):=substinpart("[",z,0), /* utility like args, but avoids / or - */)
\end{verbatim}
That's it.
\section*{EXAMPLE 1: collect terms in a sum that are constants }
\begin{verbatim}
(matchdeclare (pp, partition_expression("+",constantp,0,"+",g,'ANSpp)),
tellsimp(foo(pp),ANSpp),
declare([a,b,c],constant),
[foo(a+b+c+x+y+3),foo(w),foo(a), foo(a*b+x),foo(a*b)])
\end{verbatim}
Here the intent is to separate components of a sum into two parts:
those items that are constant versus those that are not. Each of the two
parts will be a sum, initially zero. The result will be g(yes,no), and
{\tt ANSpp} will be set to g(yes,no).
\section*{EXAMPLE 2: collect terms in a product}
Same kind of deal as the previous example, but a product.
\begin{verbatim}
(matchdeclare(qq, partition_expression("*",constantp,1,"*",g,'ANSqq)),
tellsimp(bar(qq),ANSqq),
[bar(a*b*3*x*y),bar(w),bar(a+b)])
\end{verbatim}
Here the intent is to separate components of a product into two parts:
those items that are constant versus those that are not. Each of the two
parts will be a product, initially one. The result will be g(yes,no), and
ANSqq will be set to g().
\section*{EXAMPLE 3, separate terms in a product in lists}
Same as the previous example but collecting lists of terms.
\begin{verbatim}
(matchdeclare(ss,
partition_expression("*",constantp,[],cons,"[",'ANSss)),
tellsimp(bar2(ss),ANSss),
tellsimp(bar3(ss), The_Partitions_Are(ANSss)),
[bar2(3*%pi*xx), bar2(w),bar2(%pi+x+3), bar3(x*3)])
\end{verbatim}
\section*{EXAMPLE 4, separate even and odd explicit powers of parameter \%v}
This
also illustrates the use of SETS {} for collecting data, which might
be nicer than lists for some purposes.
Note the difference between finding One odd power and finding them all.
\begin{verbatim}
(matchdeclare(odd,oddp, nov, freeof(%v)),
defmatch(OneOddPowerp,nov*%v^odd), /*%v is global */
matchdeclare(tt,
partition_expression("+",OneOddPowerp,{},adjoin,"[",'ANStt)),
tellsimp(separatepowers(tt),ANStt),
[block([%v:x],separatepowers(3*x+4*x^5+7*x^10)),
block([%v:y],separatepowers(3*y+4*y^5+7*x^10)),
block([%v:a],separatepowers(expand((a+x)^5)))])
\end{verbatim}
\section*{EXAMPLE 5: Implementing a match for f(u)du in an integral}
This program looks for a pattern that can be integrated by a simple
lookup. It accesses a list of integration formulas keyed on the name
of the function f which is discovered only after searching through the
expression. For certain functions it can compute $\int f(u)du$. In
order to do this, the program makes a list of all the factors of the
integrand, and tries some divisions.
We can use {\tt intfudu(expr,x)} directly, or if the expression is
not presented as a product, it may be worthwhile to
try {\tt intfudu(factor(expr),x)}.
This program works primarily for functions of a single argument, but
we set it up so that it can be generalized to functions of any number
of arguments. The example we give is for
exponentiation. That is, the ``\verb|^|'' operator for $u^2$ which
is encoded as {\tt expt(u,2)}, for $1/u$ which is encoded as
{\tt expt(u,-1)} or for $2^u$ which is {\tt expt(2,u)}. Also it will
therefore work for $\exp(u)$ which is {\tt expt(e,u)}.
First, set up a list of kernel functions f, which returns two items,
the integration formula and the derivative of the kernel. This latter
item is almost always the same, but for functions of several variables
like ``\verb|^|'', it matters.
\begin{verbatim}
(intable[otherwise]:=false, /*backstop for undefined kernels */
intable[sin] : lambda([u], [-cos(u),diff(u,x)]),
intable[cos] : lambda([u], [sin(u),diff(u,x)]),
intable[tan] : lambda([u], [log(sec(u)),diff(u,x)]),
intable[sec] : lambda([u], [log(tan(u)+sec(u)),diff(u,x)]),
intable[atan]: lambda([u], [-(log(u^2+1)-2*u*atan(u))/2,diff(u,x)]),
intable[log] : lambda([u], [-u+u*log(u),diff(u,x)]),
intable[sinh]: lambda([u], [cosh(u),diff(u,x)]),
intable["^"] : lambda([u,v], if freeof(x,u) then [u^v/log(u),diff(v,x)]
else if freeof(x,v)
then if (v#-1) then [u^(v+1)/(v+1), diff(u,x)]
else [log(u),diff(u,x)]),
intable[polygamma] :lambda([u,v], if freeof(x,u) then [polygamma(u-1,v),diff(v,x)]),
intable[Ci]: lambda([u], [u*Ci(u)-sin(u),diff(u,x)]),
intable[Si]: lambda([u], [u*Si(u)-cos(u),diff(u,x)]),
gradef(Ci(w), cos(w)/w),
gradef(Si(w), sin(w)/w),
intable[nounify(bessel_j)] :lambda([u,v], if(u=1) then [-bessel_j(0,v),diff(v,x)]),
intable[nounify(bessel_i)] :lambda([u,v], if(u=1) then [ bessel_i(0,v),diff(v,x)]),
intable[nounify(bessel_k)] :lambda([u,v], if(u=1) then [-bessel_k(0,v),diff(v,x)]),
/* Legendre polynomials P[n](x), presumably if n is an explicit
integer this symbolism would be removed, so we assume it is of unknown
order n. */
intable[legendre_p] :
lambda([n,u], if freeof(n,u) then
[(legendre_p(n+1,u)-legendre_(n-1,u))/(2*n+1),
diff(u,x)])
)$
/*etc etc add functions */
( matchdeclare(ss, partition_expression("*",lambda([u],freeof(x,u)),[],cons,"[",'ANSss)),
defrule(ddr1,ss,ANSss))$
intfudu(exp,x):= /*integrate exp=f(u)*du wrt x*/
block([lists,consts,factors, thefuns, therest, thelist, int, df, result:false],
if freeof(x,exp) then return (exp*x),
lists:ddr1(exp), /*partition expression into factors*/
if lists=false then lists:[[1],[exp]],
factors:second(lists),
for k in factors do
if not atom(k) and
(thefuns:intable[inpart(k,0)])#false and
(thelist:apply(thefuns,inargs(k)))
then( [int,df]:thelist,
if freeof(x,therest:ratsimp(exp/k/df))
then (result:(therest*int),
return())),
/*if nothing of form f(u)du worked, then try matching u^1*u'-> u^2/2 */
if result=false then
for k in factors do
if freeof(x,therest:ratsimp(exp/k/diff(k,x)))
then (result:(therest*k^2/2),
return()),
return(result))$
\end{verbatim}
This program appears to duplicate the first stage of Joel Moses' SIN program.
According to Moses \cite{moses}.
SIN'S first stage was able to solve 45 out of
the 86 problems presented to James Slagle's earlier SAINT program.
It probably does about 80 percent of freshman calculus problems,
excluding rational function integration,
and is likely to be rather fast. We should add some explicit examples here.
For improving this program we could add to {\tt intable}, which will not
slow the program down.
We could add more advanced stuff for functions of several variables, as shown
by the example for Bessel. Here's Polygamma (see below)
\begin{verbatim}
intable[polygamma] :lambda([u,v], if freeof(x,u) then [polygamma(u-1,v),diff(v,x)]),
\end{verbatim}
Other ways to improve could be based on Moses' stage 2, rational
function integration, and Risch integration, as well as user-contributed
special patterns, as well as heuristics such as integration by parts.
We note that while a general method may provide answers,
they may be in a form that is difficult to comprehend, so that even if ``Risch''
can do a problem, a more targeted method may help. In any case, further
work on integration goes beyond the scope of this paper.
\section*{Thanks}
Thanks to Barton Willis for running some checks and finding some bugs. All remaining bugs are, of course, the responsibility of the author.
\begin{thebibliography}{99}
\bibitem{moses}
{Joel Moses, ``Symbolic integration: the stormy decade''
Comm. ACM, Volume 14 , Issue 8 (August 1971)
548---560 }
\end{thebibliography}
\end{document}
%% some extra formulas. Cosine Integral
(intable[Ci]: lambda([u], [u*Ci[u]-sin(u),diff(u,x)]),
intable[Si]: lambda([u], [u*Si[u]-cos(u),diff(u,x)]),
gradef(Ci(w), cos(w)/w),
gradef(Si(w), sin(w)/w))
(intable[bessel_j] :lambda([u,v], if(u=1) then -bessel_j(0,v),diff(v,x)),
intable[bessel_i] :lambda([u,v], if(u=1) then bessel_i(0,v),diff(v,x)),
intable[bessel_k] :lambda([u,v], if(u=1) then -bessel_k(0,v),diff(v,x)),
/* Legendre polynomials P[n](x), presumably if n is an explicit
integer this symbolism would be removed, so we assume it is of unknown
order n. */
intable[legendre_p] :
lambda([n,u] if freeof(n,u) then
[(legendre_p(n+1,u)-legendre_(n-1,u))/(2*n+1),
diff(u,x)])
/* we could try for associated legendre polynomials, too */