I'm interested in simplifying expressions like $$\int_{-\infty}^\infty dq_1\int_{-\infty}^\infty dq_2\int_{-\infty}^\infty dq_3f(q_1,q_2,q_3)\langle\ |b(p_1)b(p_1)a^+(q_1)b(q_2)b(q_3)b^+(r_1)b^+(r_2)a^+(r_3)|\ \rangle$$ where $$a(x)|\ \rangle = b(x)|\ \rangle = 0$$ $$\langle\ |a^+(x) = \langle\ |b^+(x) = 0$$ $$[a(x),a^+(y)] = [b(x),b^+(y)] = \delta(x - y)$$ and all other commutators are 0. These arise in some quantum field theory calculations. This is doable by hand but rather tedious. I'd like to find a CAS or CAS/package combination that can do it. I've already tried the NCAlgebra package (https://mathweb.ucsd.edu/~ncalg/) for Mathematica and that works ok but I have to sort of trick it into doing what I want; it doesn't seem designed for this kind of thing. In particular, it doesn't seem to support simplification using commutators (which is typically how one simplifies expression like the one above) and Mathematica itself struggles to deal with simplifying all the Dirac delta functions that are produced.
All that being said, does any one have any recommendations for a better CAS or CAS/package combination? The ideal would be a solution where I could, for example, give it the above information and it would do the rest.
I appreciate any and all advice.