+/** An utility function looking for given metric within exvector,
+ * used in cliffordunit::contract_with() */
+int find_same_metric(exvector & v, ex & c){
+ for (int i=0; i<v.size();i++){
+ if (!is_a<clifford>(v[i]) && is_a<indexed>(v[i])
+ && ex_to<clifford>(c).same_metric(v[i])
+ && (ex_to<varidx>(c.op(1)) == ex_to<indexed>(v[i]).get_indices()[0]
+ || ex_to<varidx>(c.op(1)).toggle_variance() == ex_to<indexed>(v[i]).get_indices()[0])){
+ return ++i; // next to found
+ }
+ }
+ return 0; //nothing found
+}
+
+/** Contraction of a Clifford units with something else. */
+bool cliffordunit::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
+{
+
+ GINAC_ASSERT(is_a<clifford>(*self));
+ GINAC_ASSERT(is_a<indexed>(*other));
+ GINAC_ASSERT(is_a<cliffordunit>(self->op(0)));
+ clifford unit = ex_to<clifford>(*self);
+ unsigned char rl = unit.get_representation_label();
+
+ if (is_a<clifford>(*other)) {
+ // Contraction only makes sense if the represenation labels are equal
+ // and the metrics are the same
+ if ((ex_to<clifford>(*other).get_representation_label() != rl)
+ && unit.same_metric(*other))
+ return false;
+ // Find if a previous contraction produces the square of self
+ int prev_square = find_same_metric(v, self[0]);
+ varidx d((new symbol)->setflag(status_flags::dynallocated), ex_to<idx>(ex_to<idx>(self->op(1)).get_dim()));
+ ex squared_metric = unit.get_metric(self->op(1), d)*unit.get_metric(d.toggle_variance(), other->op(1));
+
+ // e~mu e.mu = Tr ONE
+ if (other - self == 1) {
+ if (prev_square != 0) {
+ *self = squared_metric;
+ v[prev_square-1] = _ex1;
+ } else
+ *self = unit.get_metric(self->op(1), other->op(1));
+ *other = dirac_ONE(rl);
+ return true;
+
+ // e~mu e~alpha e.mu = (2e~alpha^2-Tr) e~alpha
+ } else if (other - self == 2
+ && is_a<clifford>(self[1])) {
+
+ const ex & ia = self[1].op(1);
+ const ex & ib = self[1].op(1);
+ if (is_a<tensmetric>(unit.get_metric().op(0)))
+ *self = 2 - unit.get_metric(self->op(1), other->op(1));
+ else if (prev_square != 0) {
+ *self = 2-squared_metric;
+ v[prev_square-1] = _ex1;
+ } else
+ *self = 2*unit.get_metric(ia, ib) - unit.get_metric(self->op(1), other->op(1));
+ *other = _ex1;
+ return true;
+
+ // e~mu S e~alpha e.mu = 2 e~alpha^3 S - e~mu S e.mu e~alpha
+ // (commutate contracted indices towards each other, simplify_indexed()
+ // will re-expand and re-run the simplification)
+ } else {
+ exvector::iterator it = self + 1, next_to_last = other - 1;
+ while (it != other) {
+ if (!is_a<clifford>(*it))
+ return false;
+ ++it;
+ }
+
+ it = self + 1;
+ ex S = _ex1;
+ while (it != next_to_last) {
+ S *= *it;
+ *it++ = _ex1;
+ }
+
+ const ex & ia = next_to_last->op(1);
+ const ex & ib = next_to_last->op(1);
+ if (is_a<tensmetric>(unit.get_metric().op(0)))
+ *self = 2 * (*next_to_last) * S - (*self) * S * (*other) * (*next_to_last);
+ else if (prev_square != 0) {
+ *self = 2 * (*next_to_last) * S - (*self) * S * (*other) * (*next_to_last)*unit.get_metric(self->op(1),self->op(1));
+ v[prev_square-1] = _ex1;
+ } else
+ *self = 2 * (*next_to_last) * S* unit.get_metric(ia,ib) - (*self) * S * (*other) * (*next_to_last);
+ *next_to_last = _ex1;
+ *other = _ex1;
+ return true;
+ }
+
+ }
+
+ return false;
+}
+