sum_kron_delta -- summation of Kronecker delta in Maxima

Nov 13th, 2017
4,084
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. /* sum_kron_delta.mac -- simplify 'sum(f(i)*kron_delta(i, j), i, ...)
2.  * copyright 2013 by Robert Dodier
3.  * I release this work under terms of the GNU General Public License
4.  */
5.
6. put ('sum_kron_delta, true, 'present);
7.
8. matchdeclare (kk, kron_delta_p);
9. kron_delta_p (e) := not atom(e) and op(e) = 'kron_delta;
10. matchdeclare (xx, lambda ([e], not kron_delta_p (e)));
11. matchdeclare (ii, symbolp);
12. matchdeclare ([ii1, iin], all);
13.
14. simp : false;
15. tellsimpafter ('sum (kk, ii, ii1, iin), FOO (kk, 1, ii, ii1, iin));
16. tellsimpafter ('sum (kk * xx, ii, ii1, iin), FOO (kk, xx, ii, ii1, iin));
17. simp:true;
18.
19. FOO (kk, xx, ii, ii1, iin) := block ([ii_args],
20.   kk : kron_delta_product (kk),
21.   if atom(kk) or not (op(kk) = 'kron_delta)
22.     then 'sum (kk*xx, ii, ii1, iin)
23.     else
24.      (ii_args : sublist (args (kk), lambda ([e], not freeof (ii, e))),
25.       /* just take the first kron_delta variable other than ii -- this is not unique when there are 3 or more variables!! */
26.       block ([a],
27.         a : first (delete_all (ii_args, args (kk))),
28.         /* oops -- what if ii_args has more (or less) than 1 element?? */
29.         eq : solve (first (ii_args) = a, ii),
30.         /* oops -- what if eq has more (or less) than 1 element?? */
31.         buildq
32.         ([eq,
33.           S : index_set (ii1, iin),
34.           kk_xx : subst (eq, kk * xx),
35.           rhs_eq : rhs (first (eq))],
36.           if %elementp (rhs_eq, S) then kk_xx else 0))));
37.
38. kron_delta_product (k) := if op(k) = "*" then apply (append, args (k)) else k;
39. delete_all (L1, L2) := (for x in L1 do L2 : delete (x, L2), L2);
40. index_set (i, j) := %intersection (%range (i, j), 'integers);
41.
42. matchdeclare (xx, all);
43. matchdeclare (%ii, lambda ([e], not atom(e) and op(e) = '%intersection));
44. tellsimp (%elementp (xx, %ii), apply ("and", map (lambda ([S], %elementp (xx, S)), args (%ii))));
45.
46. matchdeclare (ii, lambda ([e], integerp (e) or featurep (e, 'integer)));
47. tellsimp (%elementp (ii, 'integers), true);
48. matchdeclare ([aa, bb], all);
49. tellsimp (%elementp (xx, %range (aa, bb)), aa <= xx and xx <= bb);