Use cyclicity of the trace

Analyzing an expression is easiest if all n! orderings of a given
product produce the same output when passed to sort_product().
Unfortunately, this is very hard to achieve in a CAS for non-commuting
objects. Consider matrices A and C with a scalar B. CAB and BCA are
equal but the first is "sorted" because B comes after A and the second
is "sorted" because B comes before C.

An opportunity for making this better occurs when the expression is
inside a trace. Then, B can be pulled out by untrace() and the rest can
be sorted with cyclicity.

Signed-off-by: Connor Behan <[email protected]>

Author: cbehan

core/algorithms/sort_product.cc

@@ -1,6 +1,7 @@
 
 #include "Cleanup.hh"
 #include "algorithms/sort_product.hh"
+#include "properties/Trace.hh"
 
 using namespace cadabra;
 
@@ -88,6 +89,61 @@ Algorithm::result_t sort_product::apply(iterator& st)
 			}
 		}
 
+	if(ret==result_t::l_no_action) {
+		bool found=false;
+		bool ascend=true;
+		iterator parn=st;
+		// It is fine to cycle when we have a sum of a sum inside a trace
+		while(ascend && !tr.is_head(parn)) {
+			parn=tr.parent(parn);
+			const Trace *trace=kernel.properties.get<Trace>(parn);
+			if(trace) {
+				ascend=false;
+				found=true;
+				}
+			else if(*parn->name=="\\indexbracket") {
+				ascend=false;
+				if(number_of_indices(parn)==2) {
+					index_iterator first=begin_index(parn);
+					index_iterator last=first;
+					++last;
+					if(*first->name==*last->name) found=true;
+					}
+				}
+			else if(*parn->name!="\\sum") {
+				ascend=false;
+				}
+			}
+		if(found) {
+			one=tr.begin(st);
+			two=one;
+			++two;
+			while(two!=tr.end(st)) {
+				compare.clear();
+				auto es=compare.equal_subtree(one, two);
+				if(compare.should_swap(one, es)) one=two;
+				++two;
+				}
+			// We have found the element that should go at the front of the trace
+			Ex::sibling_iterator front=one;
+			while(tr.begin(st)!=front) {
+				one=tr.begin(st);
+				two=one;
+				++two;
+				while(two!=tr.end(st)) {
+					compare.clear();
+					auto es=compare.equal_subtree(one, two);
+					int sign=compare.can_swap_components(one, two, es);
+					if(sign==-1) flip_sign(st->multiplier);
+					tr.swap(one);
+					++two;
+					++two;
+					}
+				}
+			ret=result_t::l_applied;
+			}
+		}
+
 	if(cleanup)
 		cleanup_dispatch(kernel, tr, st);
 

tests/implicit.cdb

@@ -144,3 +144,29 @@ def test11():
 
 test11()
 
+def test12():
+    __cdbkernel__=create_scope()
+    {a,b}::Indices(vector);
+    Tr{#}::Trace(indices=vector);
+    A::ImplicitIndex(A_{a b});
+    B::ImplicitIndex(B_{a b});
+    {A_{a b}, B_{a b}}::AntiCommuting;
+    ex:=Tr( B A );
+    sort_product(_);
+    tst:= Tr( -A B ) - @(ex);
+    assert(tst==0)
+    print("Test 12 passed")
+
+test12()
+
+def test13():
+    __cdbkernel__=create_scope()
+    {A,B,C}::ImplicitIndex;
+    ex:=(B A C)_{a a};
+    sort_product(_);
+    tst:= (A C B)_{a a} - @(ex);
+    assert(tst==0)
+    print("Test 13 passed")
+
+test13()
+