Custom simplification (rules and patterns)

Question

Hello I'm very new to maxima. How could I have something like dot products noticed and simplfied?

Attempt 1:

tellsimpafter(a[1]*b[1]+a[2]*b[2], dot(a,b));

Works for a[1]*b[1]+a[2]*b[2] but not v[1]*w[1]+v[2]*w[2] nor w[2]*v[2]+123+v[1]*w[1].

Attempt 2:

matchdeclare(a, lambda([x], length(x)=2));
matchdeclare(b, lambda([x], length(x)=2));
tellsimpafter(a[1]*b[1]+a[2]*b[2], dot(a,b));

Doesn't even work for a[1]*b[1]+a[2]*b[2].

Answer 1

Here's a possible solution. Assuming that you know the names of the arrays, it's straightforward to extract the inner product. Then loop over that with all possible array names, as determined by looking at the subscripted variables in the expression.

/* solution for https://stackoverflow.com/questions/69673365/custom-simplification-rules */

/* try all combinations of potential arrays in succession */

contract_inner_products (e) :=
 (for S in powerset (subscripted_variables(e), 2)
    do block ([%xx: first(S), %yy: second(S)],
              e: apply1 (e, rule_contract_inner_product)),
  e);

/* extract potential arrays */

subscripted_variables (e) :=
  setify (map (op, sublist (listofvars(e), subvarp)));

/* contract inner product, assuming arrays have been identified */

matchdeclare (aa, lambda([e], freeof(%xx, e) or freeof(%yy, e)),
              bb, lambda([e], freeof(%xx, e) and freeof(%yy, e)),
              xx, lambda([e], e = %xx), yy, lambda([e], e = %yy));

defrule (rule_contract_inner_product, aa + bb*xx[1]*yy[1] + bb*xx[2]*yy[2], aa + bb*dot(xx, yy));

Here are some examples.

(%i3) dotgh: g[1]*h[1] + g[2]*h[2] $
(%i4) dotab: a[1]*b[1] + a[2]*b[2] $
(%i5) dotxy: x[1]*y[1] + x[2]*y[2] $
(%i6) contract_inner_products (dotgh);
(%o6)                       dot(g, h)
(%i7) contract_inner_products (dotgh + dotxy);
(%o7)                 dot(x, y) + dot(g, h)
(%i8) contract_inner_products (dotgh - dotxy);
(%o8)                 dot(g, h) - dot(x, y)
(%i9) contract_inner_products (dotgh - 2*dotxy);
(%o9)                dot(g, h) - 2 dot(x, y)
(%i10) contract_inner_products (n/2*dotgh - 2*dotxy);
                    dot(g, h) n
(%o10)              ----------- - 2 dot(x, y)
                         2
(%i11) contract_inner_products (n/2*dotgh - 2*dotxy - 6*%pi^2);
                               dot(g, h) n        2
(%o11)       (- 2 dot(x, y)) + ----------- - 6 %pi
                                    2

It seems to work for stuff of the form dot(x, y) + dot(x, z).

(%i12) contract_inner_products (dotgh + g[1]*f[1] + g[2]*f[2]);
(%o12)                dot(g, h) + dot(f, g)
(%i13) contract_inner_products (dotgh + g[1]*f[1] + g[2]*f[2] + 123);
(%o13)             dot(g, h) + dot(f, g) + 123
(%i14) contract_inner_products (dotgh + g[1]*x[1] + g[2]*x[2] + 123);
(%o14)             dot(g, x) + dot(g, h) + 123
(%i15) contract_inner_products (dotgh + h[1]*x[1] + h[2]*x[2] + 123);
(%o15)             dot(h, x) + dot(g, h) + 123

I didn't really expect it to work for dot(x, y) + dot(x, z), but anyway that's great.

If you try it on other examples, you'll probably find some cases which this approach doesn't recognize.

EDIT: Note that it's not necessary for the would-be dot product to be at the top level of the expression. It looks for potential inner products in subexpressions, too.

(%i11) 1/(1 - dotab/4);
                                1
(%o11)                  -----------------
                            a  b  + a  b
                             2  2    1  1
                        1 - -------------
                                  4
(%i12) contract_inner_products(%);
                                1
(%o12)                    -------------
                              dot(a, b)
                          1 - ---------
                                  4
(%i13) expand (%o11);
                                1
(%o13)                ---------------------
                         a  b     a  b
                          2  2     1  1
                      (- -----) - ----- + 1
                           4        4
(%i14) contract_inner_products(%);
                                1
(%o14)                    -------------
                              dot(a, b)
                          1 - ---------
                                  4