OK, here's my attempt at a solution.
I used defrule instead of tellsimpafter to define rules.
I chose defrule in order to work around a tellsimpafter
idiosyncrasy; namely that multiple tellsimpafter rules are not all
executed -- there is a flag AFTERFLAG* or something like that
which, I guess, is supposed to prevent infinite loops, but also prevents
multiple rules for the same operator from firing.
I'd be happy to explain the rules in more detail (assuming they work
as intended). I think the built-in rule machinery is the best way to
go here -- doing stuff like factoring out scalars by hand is essentially
just reinventing stuff that already exists in the rule code.
best
Robert Dodier
PS. Here's the result of the following example. It seems plausible,
but I didn't check it. I tried some simpler examples which turned out
correct as far as I can tell.
load (commutators);
declare ([s1, s2, s3, s4], scalar);
declare ([a, b, c], nonscalar);
e : comm (3*s1*a^^2 + 4*s2*b^^3 , (s3*a^^4)/5 - 6*s4*c^^5);
FOO (e);
expand (%);
-24*s2*s4*c^^4 . comm(b,c) . b^^2-18*s1*s4*c^^4 . comm(a,c) . a
-96*s2*s4*c^^3 . comm(b,c) . c . b^^2
-72*s1*s4*c^^3 . comm(a,c) . c . a
-144*s2*s4
*c^^2 . comm(b,c) . c^^2 . b^^2
-108*s1*s4*c^^2 . comm(a,c) . c^^2 . a
-96*s2*s4*c . comm(b,c) . c^^3 . b^^2
-72*s1*s4*c . comm(a,c) . c^^3 . a
-24*s2*s4*comm(b,c) . c^^4 . b^^2
-24*s2*s4*b^^2 . c^^4 . comm(b,c)
-96*s2*s4*b^^2 . c^^3 . comm(b,c) . c
-144*s2*s4
*b^^2 . c^^2 . comm(b,c) . c^^2
-96*s2*s4*b^^2 . c . comm(b,c) . c^^3
-24*s2*s4*b^^2 . comm(b,c) . c^^4
-4*s2*s3*b^^2 . comm(a,b) . a^^3/5
-48*s2*s4*b . c^^4 . comm(b,c) . b
-192*s2*s4
*b . c^^3 . comm(b,c) . c . b
-288*s2*s4
*b . c^^2 . comm(b,c) . c^^2 . b
-192*s2*s4
*b . c . comm(b,c) . c^^3 . b
-48*s2*s4*b . comm(b,c) . c^^4 . b
-8*s2*s3*b . comm(a,b) . b . a^^3/5
-18*s1*s4*comm(a,c) . c^^4 . a
-4*s2*s3*comm(a,b) . b^^2 . a^^3/5
-4*s2*s3*a^^3 . b^^2 . comm(a,b)/5
-8*s2*s3*a^^3 . b . comm(a,b) . b/5
-4*s2*s3*a^^3 . comm(a,b) . b^^2/5
-12*s2*s3*a^^2 . b^^2 . comm(a,b) . a
/5
-24*s2*s3*a^^2 . b . comm(a,b) . b . a
/5
-12*s2*s3*a^^2 . comm(a,b) . b^^2 . a
/5-18*s1*s4*a . c^^4 . comm(a,c)
-72*s1*s4*a . c^^3 . comm(a,c) . c
-108*s1*s4*a . c^^2 . comm(a,c) . c^^2
-72*s1*s4*a . c . comm(a,c) . c^^3
-12*s2*s3*a . b^^2 . comm(a,b) . a^^2
/5
-24*s2*s3*a . b . comm(a,b) . b . a^^2
/5-18*s1*s4*a . comm(a,c) . c^^4
-12*s2*s3*a . comm(a,b) . b^^2 . a^^2
/5
PPS. Here's the code at the moment.
/* copyright 2010 by Robert Dodier.
* I release this work under terms of the GNU General Public License.
*/
/* from Ted Wollett 2010-11-20:
The ultimate objective is to symbolically evaluate things like
comm (3*s1*a^^2 + 4*s2*b^^3 , (s3*a^^4)/5 - 6*s4*c^^5),
with (s1,s2,s3) being declared scalars, and
a^^2 = a . a, for example, and
with the final result in terms of comm (a, b) = -comm (b, a),
comm (a, c) and comm (b, c), letting the user then
replace the basic commutators with whatever the algebra
dictates.
The use of multiadditive gets the linear in each arg property
of the commutator.
In addition to the antisymmetric property of comm, we need
to obey the basic product rules:
comm (a . b ,c) --> a . comm(b, c) + comm (a, c) . b
comm (a, b . c) --> b . comm (a, c) + comm (a, b) . c
in which the dot "." is being used for noncommutative
multiplication, since (a, b, c) here are symbols which
represent operators, matrices, etc., and hence we need
to preserve the correct ordering of (a,b,c) in the
final result.
*/
declare (bilinear, feature);
bilinearp (e) := featurep (e, bilinear);
declare (comm, bilinear);
matchdeclare (xx, bilinearp);
matchdeclare ([aa, bb, cc], all);
defrule (rbilineara1, xx (aa + bb, cc), map (lambda ([u], xx (u, cc)), bb));
defrule (rbilineara2, xx (cc, aa + bb), map (lambda ([u], xx (cc, u)), bb));
matchdeclare (xx, bilinearp);
matchdeclare (aa, scalarp);
matchdeclare (bb, nonscalarp);
matchdeclare (cc, all);
defrule (rbilinearm1, xx (aa * bb, cc), aa * xx (bb, cc));
defrule (rbilinearm2, xx (cc, aa * bb), aa * xx (cc, bb));
nctimesp (e) := not atom(e) and op(e) = ".";
matchdeclare (aa, nctimesp);
matchdeclare (cc, all);
defrule (rcommp1, comm (aa, cc), first (aa) . comm (rest (aa), cc) +
comm (rest (aa, -1), cc) . last (aa));
defrule (rcommp2, comm (cc, aa), first (aa) . comm (cc, rest (aa)) +
comm (cc, rest (aa, -1)) . last (aa));
matchdeclare (aa, nonscalarp);
matchdeclare (nn, integerp);
defrule (rncexpt, aa^^nn, apply (".", makelist (aa, i, 1, nn)));
/* put antisymmetric declaration here so built-in rules don't
interfere w/ rule definitions */
declare (comm, antisymmetric);
FOO (e) :=
(e : apply1 (e, rbilineara1, rbilineara2, rbilinearm1, rbilinearm2),
block ([dotexptsimp : false],
apply1 (applyb1 (e, rncexpt), rcommp1, rcommp2)));
/* example:
*
load (commutators);
declare ([s1, s2, s3, s4], scalar);
declare ([a, b, c], nonscalar);
e : comm (3*s1*a^^2 + 4*s2*b^^3 , (s3*a^^4)/5 - 6*s4*c^^5);
FOO (e);
expand (%);
*
*/